(*     Solving Large(n > 2), Linear Systems of n Equations in n Hypercomplex Unknowns - Using The Eight-Step Method

                                      by Chris Niemann  *)

(*  The solution is based on two key ingredients: the Transformation Matrix and the Eight-Step Method. The two 
Transformation Matrices are formed by computing the product from the left and the product from the
right. Both of these matrices are needed because one could have a constant on the left and/or the right of the unknown. First we illustrate the fact that the same comstants on the left versus on the right will lead to two distinct solutions.  *)

QuatDot[x_,y_]:=(*  User must "ENTER" here.  *)      (*  This defines the 4-dimensional hypercomplex product.  *)(If[Length[x]≠4||Length[y]≠4,Return["Undefined"] ,      Return[{x[[1]]y[[1]]-x[[2]]y[[2]]-x[[3]]y[[3]]-x[[4]]y[[4]],x[[1]]y[[2]]+x[[2]]y[[1]]+x[[3]]y[[4]]-x[[4]]y[[3]],x[[1]]y[[3]]-x[[2]]y[[4]]+x[[3]]y[[1]]+x[[4]]y[[2]],x[[1]]y[[4]]+x[[2]]y[[3]]-x[[3]]y[[2]]+x[[4]]y[[1]]}]]);QuatSqr[x_]:=(QuatDot[x,x]);

FirstHalf[q_]:=(Lq=Length[q];L2=Lq/2;Tyler={q[[1]]};Do[Tyler=Append[Tyler,q[[k]]],{k,2,L2}];Return[Tyler]);
SecondHalf[q_]:=(Lq=Length[q];L2=Lq/2;Kaleb={q[[L2+1]]};Do[Kaleb=Append[Kaleb,q[[k]]],{k,L2+2,Lq}];Return[Kaleb]);

Hconj[q_]:=(tmp=-q;tmp[[1]]=-tmp[[1]];Return[tmp]);  (*  Returns the hypercomplex conjugate of  q.  *)

(*  Below makes use of the Caley-Dickson product fromula to compute the octonion product of  x  followed by  y.  *)
OctDot[x_,y_]:=(If[Length[x]≠8||Length[y]≠8,Return["Undefined"]];(Subscript[Amanda, 4])=FirstHalf[x];(Subscript[Hannah, 4])=SecondHalf[x];
(Subscript[Zoe, 4])=FirstHalf[y];(Subscript[Ainsley, 4])=SecondHalf[y];(Subscript[Kyle, 4])=QuatDot[(Subscript[Amanda, 4]),(Subscript[Zoe, 4])]-QuatDot[Hconj[(Subscript[Ainsley, 4])],(Subscript[Hannah, 4])];(Subscript[Cole, 4])=QuatDot[(Subscript[Hannah, 4]),Hconj[(Subscript[Zoe, 4])]]+QuatDot[(Subscript[Ainsley, 4]),(Subscript[Amanda, 4])];Return[Join[(Subscript[Kyle, 4]),(Subscript[Cole, 4])]]);OctSqr[x_]:=(OctDot[x,x]);

(*  Below makes use of the Caley-Dickson product fromula to compute the sedenion product of  x  followed by  y.  *)SedDot[x_,y_]:=(If[Length[x]≠16||Length[y]≠16,Return["Undefined"]];(Subscript[Amanda, 8])=FirstHalf[x];(Subscript[Hannah, 8])=SecondHalf[x];
(Subscript[Zoe, 8])=FirstHalf[y];(Subscript[Ainsley, 8])=SecondHalf[y];(Subscript[Kyle, 8])=OctDot[(Subscript[Amanda, 8]),(Subscript[Zoe, 8])]-OctDot[Hconj[(Subscript[Ainsley, 8])],(Subscript[Hannah, 8])];(Subscript[Cole, 8])=OctDot[(Subscript[Hannah, 8]),Hconj[(Subscript[Zoe, 8])]]+OctDot[(Subscript[Ainsley, 8]),(Subscript[Amanda, 8])];Return[Join[(Subscript[Kyle, 8]),(Subscript[Cole, 8])]]);SedSqr[x_]:=(SedDot[x,x])
Quarter4x4[M_]:=((Subscript[Q, 1])={{M[[1,1]],M[[1,2]]},{M[[2,1]],M[[2,2]]}};(Subscript[Q, 2])={{M[[1,3]],M[[1,4]]},{M[[2,3]],M[[2,4]]}};(Subscript[Q, 3])={{M[[3,1]],M[[3,2]]},{M[[4,1]],M[[4,2]]}};(Subscript[Q, 4])={{M[[3,3]],M[[3,4]]},{M[[4,3]],M[[4,4]]}});

(*   To quote from the wikipedia web site,"In mathematics,the general linear group of degree n is the set of n×n invertible matrices,together with the operation of ordinary matrix multiplication. This forms a group because the product of two invertible matrices is again invertible and the inverse of an invertible matrix is invertible. The name is because the columns of an invertible matrix are linearly independent, hence the vectors/points they define are in general linear position, and matrices in the general linear group take points in general linear position to points in general linear position."  The General Linear Group of degree n over any field, F ->  (GL(n,F)  is defined below  *)
GL[n_,F_]:=(If[F≤1||n≤0,Return["Undefined"]];GLG=1;Fpwr=0;Label[continue];factor=F^n-F^Fpwr;GLG=GLG*factor;Fpwr=Fpwr+1;If[Fpwr==n,Return[GLG],Goto[continue]]);

(*   Below are the four Eight-Step Methods developed by Chirs Niemann from 1988-1990+  Any of the four may be used as long as the matrix of that quarter ((Subscript[Q, 1]) -> A, (Subscript[Q, 2]) -> B, (Subscript[Q, 3]) -> C & (Subscript[Q, 4]) -> D) is invertible. The matrix ({{A, B}, {C, D}}) is invertible iff the matrix (Subscript[T, n]) (for n = 1, 2, 3, or 4) is nonsingular. (In PEP4, T in known as Schur's complement). These methods are particularly useful in solving linear systems of hypercpmplex numbers because the four quarters are naturally formed via the transformation matrices.     *) (Subscript[PEP, 1])[A_,B_,C_,D_]:=((Subscript[R, 1])=-Inverse[A];(Subscript[S, 1])=(Subscript[R, 1]).B;(Subscript[T, 1])=C.(Subscript[S, 1])+D;Z=Inverse[(Subscript[T, 1])];X=(Subscript[S, 1]).Z;(Subscript[V, 1])=C.(Subscript[R, 1]);Y=Z.(Subscript[V, 1]);W=(Subscript[S, 1]).Y-(Subscript[R, 1]);Return[]);(Subscript[PEP, 2])[A_,B_,C_,D_]:=((Subscript[R, 2])=-Inverse[B];(Subscript[S, 2])=(Subscript[R, 2]).A;(Subscript[T, 2])=D.(Subscript[S, 2])+C;X=Inverse[(Subscript[T, 2])];Z=(Subscript[S, 2]).X;(Subscript[V, 2])=D.(Subscript[R, 2]);W=X.(Subscript[V, 2]);Y=(Subscript[S, 2]).W-(Subscript[R, 2]);Return[]);(Subscript[PEP, 3])[A_,B_,C_,D_]:=((Subscript[R, 3])=-Inverse[C];(Subscript[S, 3])=(Subscript[R, 3]).D;(Subscript[T, 3])=A.(Subscript[S, 3])+B;Y=Inverse[(Subscript[T, 3])];W=(Subscript[S, 3]).Y;(Subscript[V, 3])=A.(Subscript[R, 3]);Z=Y.(Subscript[V, 3]);X=(Subscript[S, 3]).Z-(Subscript[R, 3]);Return[]);(Subscript[PEP, 4])[A_,B_,C_,D_]:=((Subscript[R, 4])=-Inverse[D];(Subscript[S, 4])=(Subscript[R, 4]).C;(Subscript[T, 4])=B.(Subscript[S, 4])+A;W=Inverse[(Subscript[T, 4])];Y=(Subscript[S, 4]).W;(Subscript[V, 4])=B.(Subscript[R, 4]);X=W.(Subscript[V, 4]);Z=(Subscript[S, 4]).X-(Subscript[R, 4]);Return[]);

(Subscript[TransMatL, 4])[q_]:=  (*  This is the Transformation Matrix corresponding to quaternion multiplication from the Left.   *)({{q[[1]],-q[[2]],-q[[3]],-q[[4]]},{q[[2]],q[[1]],-q[[4]],q[[3]]},{q[[3]],q[[4]],q[[1]],-q[[2]]},{q[[4]],-q[[3]],q[[2]],q[[1]]}});(Subscript[TransMatR, 4])[q_]:= (*  The Transformation Matrix below corresponds to quaternion multiplication from the Right.   *)({{q[[1]],-q[[2]],-q[[3]],-q[[4]]},{q[[2]],q[[1]],q[[4]],-q[[3]]},{q[[3]],-q[[4]],q[[1]],q[[2]]},{q[[4]],q[[3]],-q[[2]],q[[1]]}});(Subscript[TransMatL, 8])[q_]:=  (*  This is the Transformation Matrix corresponds to qz, where z is the octonion unknown   *)({{q[[1]],-q[[2]],-q[[3]],-q[[4]],-q[[5]],-q[[6]],-q[[7]],-q[[8]]},{q[[2]],q[[1]],-q[[4]],q[[3]],-q[[6]],q[[5]],q[[8]],-q[[7]]},{q[[3]],q[[4]],q[[1]],-q[[2]],-q[[7]],-q[[8]],q[[5]],q[[6]]},{q[[4]],-q[[3]],q[[2]],q[[1]],-q[[8]],q[[7]],-q[[6]],q[[5]]},
{q[[5]],q[[6]],q[[7]],q[[8]],q[[1]],-q[[2]],-q[[3]],-q[[4]]},
{q[[6]],-q[[5]],q[[8]],-q[[7]],q[[2]],q[[1]],q[[4]],-q[[3]]},
{q[[7]],-q[[8]],-q[[5]],q[[6]],q[[3]],-q[[4]],q[[1]],q[[2]]},
{q[[8]],q[[7]],-q[[6]],-q[[5]],q[[4]],q[[3]],-q[[2]],q[[1]]}});(Subscript[TransMatR, 8])[q_]:= (*  This is the Transformation Matrix corresponds to zq, where z is the octonion unknown   *)({{q[[1]],-q[[2]],-q[[3]],-q[[4]],-q[[5]],-q[[6]],-q[[7]],-q[[8]]},{q[[2]],q[[1]],q[[4]],-q[[3]],q[[6]],-q[[5]],-q[[8]],q[[7]]},
{q[[3]],-q[[4]],q[[1]],q[[2]],q[[7]],q[[8]],-q[[5]],-q[[6]]},
{q[[4]],q[[3]],-q[[2]],q[[1]],q[[8]],-q[[7]],q[[6]],-q[[5]]},
{q[[5]],-q[[6]],-q[[7]],-q[[8]],q[[1]],q[[2]],q[[3]],q[[4]]},
{q[[6]],q[[5]],-q[[8]],q[[7]],-q[[2]],q[[1]],-q[[4]],q[[3]]},
{q[[7]],q[[8]],q[[5]],-q[[6]],-q[[3]],q[[4]],q[[1]],-q[[2]]},
{q[[8]],-q[[7]],q[[6]],q[[5]],-q[[4]],-q[[3]],q[[2]],q[[1]]}});(Subscript[TransMatL, 16])[q_]:=  (*  This is the Transformation Matrix corresponding to sedenion multiplication from the Left.   *)({{q[[01]],-q[[02]],-q[[03]],-q[[04]],-q[[05]],-q[[06]],-q[[07]],-q[[08]],-q[[09]],-q[[10]],-q[[11]],-q[[12]],-q[[13]],-q[[14]],-q[[15]],-q[[16]]},{q[[02]],+q[[01]],-q[[04]],+q[[03]],-q[[06]],+q[[05]],+q[[08]],-q[[07]],-q[[10]],+q[[09]],+q[[12]],-q[[11]],+q[[14]],-q[[13]],-q[[16]],+q[[15]]},                                                                                    {q[[03]],+q[[04]],+q[[01]],-q[[02]],-q[[07]],-q[[08]],+q[[05]],+q[[06]],-q[[11]],-q[[12]],+q[[09]],+q[[10]],+q[[15]],+q[[16]],-q[[13]],-q[[14]]},{q[[04]],-q[[03]],+q[[02]],+q[[01]],-q[[08]],+q[[07]],-q[[06]],+q[[05]],-q[[12]],+q[[11]],-q[[10]],+q[[09]],+q[[16]],-q[[15]],+q[[14]],-q[[13]]},{q[[05]],+q[[06]],+q[[07]],+q[[08]],+q[[01]],-q[[02]],-q[[03]],-q[[04]],-q[[13]],-q[[14]],-q[[15]],-q[[16]],+q[[09]],+q[[10]],+q[[11]],+q[[12]]},{q[[06]],-q[[05]],+q[[08]],-q[[07]],+q[[02]],+q[[01]],+q[[04]],-q[[03]],-q[[14]],+q[[13]],-q[[16]],+q[[15]],-q[[10]],+q[[09]],-q[[12]],+q[[11]]},{q[[07]],-q[[08]],-q[[05]],+q[[06]],+q[[03]],-q[[04]],+q[[01]],+q[[02]],-q[[15]],+q[[16]],+q[[13]],-q[[14]],-q[[11]],+q[[12]],+q[[09]],-q[[10]]},{q[[08]],+q[[07]],-q[[06]],-q[[05]],+q[[04]],+q[[03]],-q[[02]],+q[[01]],-q[[16]],-q[[15]],+q[[14]],+q[[13]],-q[[12]],-q[[11]],+q[[10]],+q[[09]]},{q[[09]],+q[[10]],+q[[11]],+q[[12]],+q[[13]],+q[[14]],+q[[15]],+q[[16]],+q[[01]],-q[[02]],-q[[03]],-q[[04]],-q[[05]],-q[[06]],-q[[07]],-q[[08]]},{q[[10]],-q[[09]],+q[[12]],-q[[11]],+q[[14]],-q[[13]],-q[[16]],+q[[15]],+q[[02]],+q[[01]],+q[[04]],-q[[03]],+q[[06]],-q[[05]],-q[[08]],+q[[07]]},{q[[11]],-q[[12]],-q[[09]],+q[[10]],+q[[15]],+q[[16]],-q[[13]],-q[[14]],+q[[03]],-q[[04]],+q[[01]],+q[[02]],+q[[07]],+q[[08]],-q[[05]],-q[[06]]},{q[[12]],+q[[11]],-q[[10]],-q[[09]],+q[[16]],-q[[15]],+q[[14]],-q[[13]],+q[[04]],+q[[03]],-q[[02]],+q[[01]],+q[[08]],-q[[07]],+q[[06]],-q[[05]]},{q[[13]],-q[[14]],-q[[15]],-q[[16]],-q[[09]],+q[[10]],+q[[11]],+q[[12]],+q[[05]],-q[[06]],-q[[07]],-q[[08]],+q[[01]],+q[[02]],+q[[03]],+q[[04]]},{q[[14]],+q[[13]],-q[[16]],+q[[15]],-q[[10]],-q[[09]],-q[[12]],+q[[11]],+q[[06]],+q[[05]],-q[[08]],+q[[07]],-q[[02]],+q[[01]],-q[[04]],+q[[03]]},{q[[15]],+q[[16]],+q[[13]],-q[[14]],-q[[11]],+q[[12]],-q[[09]],-q[[10]],+q[[07]],+q[[08]],+q[[05]],-q[[06]],-q[[03]],+q[[04]],+q[[01]],-q[[02]]},{q[[16]],-q[[15]],+q[[14]],+q[[13]],-q[[12]],-q[[11]],+q[[10]],-q[[09]],+q[[08]],-q[[07]],+q[[06]],+q[[05]],-q[[04]],-q[[03]],+q[[02]],+q[[01]]}});(Subscript[TransMatR, 16])[q_]:= (*  This is the Transformation Matrix corresponding to sedenion multiplication from the Right.   *)  ({{q[[01]],-q[[02]],-q[[03]],-q[[04]],-q[[05]],-q[[06]],-q[[07]],-q[[08]],-q[[09]],-q[[10]],-q[[11]],-q[[12]],-q[[13]],-q[[14]],-q[[15]],-q[[16]]},{q[[02]],+q[[01]],+q[[04]],-q[[03]],+q[[06]],-q[[05]],-q[[08]],+q[[07]],+q[[10]],-q[[09]],-q[[12]],+q[[11]],-q[[14]],+q[[13]],+q[[16]],-q[[15]]},{q[[03]],-q[[04]],+q[[01]],+q[[02]],+q[[07]],+q[[08]],-q[[05]],-q[[06]],+q[[11]],+q[[12]],-q[[09]],-q[[10]],-q[[15]],-q[[16]],+q[[13]],+q[[14]]},{q[[04]],+q[[03]],-q[[02]],+q[[01]],+q[[08]],-q[[07]],+q[[06]],-q[[05]],+q[[12]],-q[[11]],+q[[10]],-q[[09]],-q[[16]],+q[[15]],-q[[14]],+q[[13]]},{q[[05]],-q[[06]],-q[[07]],-q[[08]],+q[[01]],+q[[02]],+q[[03]],+q[[04]],+q[[13]],+q[[14]],+q[[15]],+q[[16]],-q[[09]],-q[[10]],-q[[11]],-q[[12]]},{q[[06]],+q[[05]],-q[[08]],+q[[07]],-q[[02]],+q[[01]],-q[[04]],+q[[03]],+q[[14]],-q[[13]],+q[[16]],-q[[15]],+q[[10]],-q[[09]],+q[[12]],-q[[11]]},{q[[07]],+q[[08]],+q[[05]],-q[[06]],-q[[03]],+q[[04]],+q[[01]],-q[[02]],+q[[15]],-q[[16]],-q[[13]],+q[[14]],+q[[11]],-q[[12]],-q[[09]],+q[[10]]},{q[[08]],-q[[07]],+q[[06]],+q[[05]],-q[[04]],-q[[03]],+q[[02]],+q[[01]],+q[[16]],+q[[15]],-q[[14]],-q[[13]],+q[[12]],+q[[11]],-q[[10]],-q[[09]]},{q[[09]],-q[[10]],-q[[11]],-q[[12]],-q[[13]],-q[[14]],-q[[15]],-q[[16]],+q[[01]],+q[[02]],+q[[03]],+q[[04]],+q[[05]],+q[[06]],+q[[07]],+q[[08]]},{q[[10]],+q[[09]],-q[[12]],+q[[11]],-q[[14]],+q[[13]],+q[[16]],-q[[15]],-q[[02]],+q[[01]],-q[[04]],+q[[03]],-q[[06]],+q[[05]],+q[[08]],-q[[07]]},{q[[11]],+q[[12]],+q[[09]],-q[[10]],-q[[15]],-q[[16]],+q[[13]],+q[[14]],-q[[03]],+q[[04]],+q[[01]],-q[[02]],-q[[07]],-q[[08]],+q[[05]],+q[[06]]},{q[[12]],-q[[11]],+q[[10]],+q[[09]],-q[[16]],+q[[15]],-q[[14]],+q[[13]],-q[[04]],-q[[03]],+q[[02]],+q[[01]],-q[[08]],+q[[07]],-q[[06]],+q[[05]]},{q[[13]],+q[[14]],+q[[15]],+q[[16]],+q[[09]],-q[[10]],-q[[11]],-q[[12]],-q[[05]],+q[[06]],+q[[07]],+q[[08]],+q[[01]],-q[[02]],-q[[03]],-q[[04]]},{q[[14]],-q[[13]],+q[[16]],-q[[15]],+q[[10]],+q[[09]],+q[[12]],-q[[11]],-q[[06]],-q[[05]],+q[[08]],-q[[07]],+q[[02]],+q[[01]],+q[[04]],-q[[03]]},{q[[15]],-q[[16]],-q[[13]],+q[[14]],+q[[11]],-q[[12]],+q[[09]],+q[[10]],-q[[07]],-q[[08]],-q[[05]],+q[[06]],+q[[03]],-q[[04]],+q[[01]],+q[[02]]},{q[[16]],q[[15]],-q[[14]],-q[[13]],+q[[12]],+q[[11]],-q[[10]],+q[[09]],-q[[08]],+q[[07]],-q[[06]],-q[[05]],+q[[04]],+q[[03]],-q[[02]],+q[[01]]}});BuildUp[n_,A_,B_,C_,D_]:=(*  This patches together 4 (n x n) quarters to return a (2n x 2n) matrix.  *)(Return[Join[Table[Join[Table[A[[i,j]],{j,n}],Table[B[[i,j]],{j,n}]],{i,n}],Table[Join[Table[C[[i,j]],{j,n}],Table[D[[i,j]],{j,n}]],{i,n}]]]);  (*  May I Join your Table?  *)
TFR[d_]:=(Table[Floor[19Random[]-9],{d}]);ZeroMat[n_]:=(Return[IdentityMatrix[n]-IdentityMatrix[n]]);(Subscript[q, 0])={0,0,0,0};(Subscript[e, 0])={1,0,0,0};

{GL[2,2],GL[3,2],GL[4,2]}(*  The number of symmetries in the fano line, the fano plane & the fano tetrahedron 
are computed here.  *)

(*  All 16 constants and the 2 unknowns below are quaternions. We solve by transforming the 2 equations in 2 quaternion unknowns into 8 equations in 8 real unknowns. The equations are given in their most general form, that is, multiplication from the left and multiplication from the right. If a particular system has only multiplication from the left, we need only use the transformation matrix "TransMat4L." Likewise if there is only a constant on the right, use "TransMat4R." Of course one could also assign a "missing" constant the 
quaternion value: 
							(Subscript[e, 0]) = {1,0,0,0}, that is, "1."
  *)

(*  To illustrate the solution we randomly select 16 quaternion constants with +/- digit components .  *)

SeedRandom[9797];(Subscript[a, 1])=TFR[4];(Subscript[a, 2])=TFR[4];TFR[4];(Subscript[b, 1])=TFR[4];(Subscript[b, 2])=TFR[4];(Subscript[c, 1])=TFR[4];(Subscript[c, 2])=TFR[4];(Subscript[d, 1])=TFR[4];(Subscript[d, 2])=TFR[4];(Subscript[a, 3])=TFR[4];(Subscript[a, 4])=TFR[4];(Subscript[b, 3])=TFR[4];(Subscript[b, 4])=TFR[4];(Subscript[c, 3])=TFR[4];(Subscript[c, 4])=TFR[4];(Subscript[d, 3])=TFR[4];(Subscript[d, 4])=TFR[4];
(*  User must "ENTER" here. You may (of course) change the seed number above to generate new systems.  *)

(*   As is often done in textbook examples, randomly select the answers so that the solution will avoid 
those ugly fractions.  *)

(*  "Answers" for (Subscript[x, 1]), ..., (Subscript[x, 4])  *)(*  User must "ENTER" below.  *)

(Subscript[x, 1])=TFR[4] ;(Subscript[x, 2])=TFR[4] ;(Subscript[x, 3])=TFR[4] ;(Subscript[x, 4])=TFR[4];

(*  Compute the rhs constants for multiplication from the left.  *)

(Subscript[k, 1])=QuatDot[(Subscript[a, 1]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 1]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 1]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 1]),(Subscript[x, 4])]           (*  User must "ENTER" here.  *)

(Subscript[k, 2])=QuatDot[(Subscript[a, 2]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 2]),(Subscript[x, 2])] +QuatDot[(Subscript[c, 2]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 2]),(Subscript[x, 4])]           (*  User must "ENTER" here.  *)

(Subscript[k, 3])=QuatDot[(Subscript[a, 3]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 3]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 3]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 3]),(Subscript[x, 4])]            (*  User must "ENTER" here.  *)

(Subscript[k, 4])=QuatDot[(Subscript[a, 4]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 4]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 4]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 4]),(Subscript[x, 4])]            (*  User must "ENTER" here.  *)

(*     We solve the quaternion system:
	 (Subscript[a, 1])(Subscript[x, 1]) + (Subscript[b, 1])(Subscript[x, 2]) + (Subscript[c, 1])(Subscript[x, 3]) + (Subscript[d, 1])(Subscript[x, 4]) = (Subscript[k, 1]),
		(Subscript[a, 2])(Subscript[x, 1]) + (Subscript[b, 2])(Subscript[x, 2]) + (Subscript[c, 2])(Subscript[x, 3]) + (Subscript[d, 2])(Subscript[x, 4]) = (Subscript[k, 2]),
		(Subscript[a, 3])(Subscript[x, 1]) + (Subscript[b, 3])(Subscript[x, 2]) + (Subscript[c, 3])(Subscript[x, 3]) + (Subscript[d, 3])(Subscript[x, 4]) = (Subscript[k, 1]),
		(Subscript[a, 4])(Subscript[x, 1]) + (Subscript[b, 4])(Subscript[x, 2]) + (Subscript[c, 4])(Subscript[x, 3]) + (Subscript[d, 4])(Subscript[x, 4]) = (Subscript[k, 1]),		
  *)

A4=(Subscript[TransMatL, 4])[(Subscript[a, 1])];B4=(Subscript[TransMatL, 4])[(Subscript[b, 1])];C4=(Subscript[TransMatL, 4])[(Subscript[a, 2])];D4=(Subscript[TransMatL, 4])[(Subscript[b, 2])];AA=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[c, 1])];B4=(Subscript[TransMatL, 4])[(Subscript[d, 1])];C4=(Subscript[TransMatL, 4])[(Subscript[c, 2])];D4=(Subscript[TransMatL, 4])[(Subscript[d, 2])];BB=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[a, 3])];B4=(Subscript[TransMatL, 4])[(Subscript[b, 3])];C4=(Subscript[TransMatL, 4])[(Subscript[a, 4])];D4=(Subscript[TransMatL, 4])[(Subscript[b, 4])];CC=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[c, 3])];B4=(Subscript[TransMatL, 4])[(Subscript[d, 3])];C4=(Subscript[TransMatL, 4])[(Subscript[c, 4])];D4=(Subscript[TransMatL, 4])[(Subscript[d, 4])];DD=BuildUp[4,A4,B4,C4,D4]

(Subscript[PEP, 1])[AA,BB,CC,DD]

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[8]&&AA.X+BB.Z==ZeroMat[8]&&CC.W+DD.Y==ZeroMat[8]&&CC.X+DD.Z==IdentityMatrix[8]

(Subscript[x, 1])

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]

(Subscript[x, 2])

Y.Join[(Subscript[k, 1]),(Subscript[k, 2])]+Z.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 3]),(Subscript[x, 4])]

(*  Of course if one enounters a 3x3 system, then merely let 
(Subscript[a, 4]) = (Subscript[b, 4]) = (Subscript[c, 4]) = (Subscript[d, 1]) = (Subscript[d, 2]) = (Subscript[d, 3]) = (Subscript[q, 0])  &  (Subscript[d, 4]) = (Subscript[e, 0]), then  (Subscript[x, 4]) = (Subscript[k, 4 ])& let (Subscript[k, 4]) = (Subscript[q, 0]). Example:
     *)

(Subscript[a, 4])=(Subscript[q, 0]);(Subscript[b, 4])=(Subscript[q, 0]);(Subscript[c, 4])=(Subscript[q, 0]);(Subscript[d, 1])=(Subscript[q, 0]);(Subscript[d, 2])=(Subscript[q, 0]);(Subscript[d, 3])=(Subscript[q, 0]);(Subscript[k, 4])=(Subscript[q, 0]);(Subscript[d, 4]) = (Subscript[e, 0]);(Subscript[x, 4]) = (Subscript[k, 4 ])

(Subscript[k, 1])=QuatDot[(Subscript[a, 1]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 1]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 1]),(Subscript[x, 3])]         (*  User must "ENTER" here.  *)

(Subscript[k, 2])=QuatDot[(Subscript[a, 2]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 2]),(Subscript[x, 2])] +QuatDot[(Subscript[c, 2]),(Subscript[x, 3])]           (*  User must "ENTER" here.  *)

(Subscript[k, 3])=QuatDot[(Subscript[a, 3]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 3]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 3]),(Subscript[x, 3])]            (*  User must "ENTER" here.  *)

A4      (*  remains unchanged  *)

A4=(Subscript[TransMatL, 4])[(Subscript[c, 1])];B4=(Subscript[TransMatL, 4])[(Subscript[d, 1])];C4=(Subscript[TransMatL, 4])[(Subscript[c, 2])];D4=(Subscript[TransMatL, 4])[(Subscript[d, 2])];BB=BuildUp[4,A4,B4,C4,D4]

MatrixForm[BB]

A4=(Subscript[TransMatL, 4])[(Subscript[a, 3])];B4=(Subscript[TransMatL, 4])[(Subscript[b, 3])];C4=(Subscript[TransMatL, 4])[(Subscript[a, 4])];D4=(Subscript[TransMatL, 4])[(Subscript[b, 4])];CC=BuildUp[4,A4,B4,C4,D4]

MatrixForm[CC]

A4=(Subscript[TransMatL, 4])[(Subscript[c, 3])];B4=(Subscript[TransMatL, 4])[(Subscript[d, 3])];C4=(Subscript[TransMatL, 4])[(Subscript[c, 4])];D4=(Subscript[TransMatL, 4])[(Subscript[d, 4])];DD=BuildUp[4,A4,B4,C4,D4]

Inverse[%]

MatrixForm[DD]

(Subscript[PEP, 1])[AA,BB,CC,DD]

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 1]),(Subscript[x, 2])]

Y.Join[(Subscript[k, 1]),(Subscript[k, 2])]+Z.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 3]),(Subscript[x, 4])]

(*  With this setup, neither (Subscript[PEP, 2]), nor (Subscript[PEP, 3]) could be used, because BB and CC are singular.  *)

(Subscript[PEP, 4])[AA,BB,CC,DD]

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 1]),(Subscript[x, 2])]&&Y.Join[(Subscript[k, 1]),(Subscript[k, 2])]+Z.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 3]),(Subscript[x, 4])]

(*  In case one insists on using (Subscript[PEP, 3]), then do an Infinitesimal Perturbation of the 3'rd quarter, CC, by inserting an
infinitesimal, like say, ϵ. 
  *)

Do[CC[[i,i]]=ϵ,{i,5,8}]   (*  We insert ϵ's into the lower right quarter's diagonal to make CC invertible.  *)

MatrixForm[CC]

(Subscript[PEP, 3])[AA,BB,CC,DD]

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 1]),(Subscript[x, 2])]&&(Y.Join[(Subscript[k, 1]),(Subscript[k, 2])]+Z.Join[(Subscript[k, 3]),(Subscript[k, 4])]/.ϵ->0)==Join[(Subscript[x, 3]),(Subscript[x, 4])]

(*  Example of using ϵ to perform diff., where  f(x) = (x^4).  *)   Simplify[((Expand[(((x+ϵ))^4)]-(x^4))/ϵ)]

%/. ϵ->0       (*  Smuggle in 0 (it's for real) in place of ϵ.  *)

(*  Of course, if BB is made invertible, then  (Subscript[PEP, 2])\IP  may be performed, as well.   *)

(*  -----------------------------------  fin: 4x4 Quaternion matrices ------------------------------------------------  *)

(Subscript[a, 1])=TFR[4];(Subscript[a, 2])=TFR[4];(Subscript[b, 1])=TFR[4];(Subscript[b, 2])=TFR[4];(Subscript[c, 1])=TFR[4];(Subscript[c, 2])=TFR[4];(Subscript[d, 1])=TFR[4];(Subscript[d, 2])=TFR[4];(Subscript[a, 3])=TFR[4];(Subscript[a, 4])=TFR[4];(Subscript[b, 3])=TFR[4];(Subscript[b, 4])=TFR[4];(Subscript[c, 3])=TFR[4];(Subscript[c, 4])=TFR[4];(Subscript[d, 3])=TFR[4];(Subscript[d, 4])=TFR[4];(Subscript[a, 5])=TFR[4];(Subscript[a, 6])=TFR[4];(Subscript[b, 5])=TFR[4];(Subscript[b, 6])=TFR[4];(Subscript[c, 5])=TFR[4];(Subscript[c, 6])=TFR[4];(Subscript[d, 5])=TFR[4];(Subscript[d, 6])=TFR[4];(Subscript[a, 7])=TFR[4];(Subscript[a, 8])=TFR[4];(Subscript[b, 7])=TFR[4];(Subscript[b, 8])=TFR[4];(Subscript[c, 7])=TFR[4];(Subscript[c, 8])=TFR[4];(Subscript[d, 7])=TFR[4];(Subscript[d, 8])=TFR[4];

(Subscript[e, 1])=TFR[4];(Subscript[f, 1])=TFR[4];(Subscript[g, 1])=TFR[4];(Subscript[h, 1])=TFR[4];(Subscript[e, 2])=TFR[4];(Subscript[f, 2])=TFR[4];(Subscript[g, 2])=TFR[4];(Subscript[h, 2])=TFR[4];(Subscript[e, 3])=TFR[4];(Subscript[f, 3])=TFR[4];(Subscript[g, 3])=TFR[4];(Subscript[h, 3])=TFR[4];(Subscript[e, 4])=TFR[4];(Subscript[f, 4])=TFR[4];(Subscript[g, 4])=TFR[4];(Subscript[h, 4])=TFR[4];(Subscript[e, 5])=TFR[4];(Subscript[f, 5])=TFR[4];(Subscript[g, 5])=TFR[4];(Subscript[h, 5])=TFR[4];(Subscript[e, 6])=TFR[4];(Subscript[f, 6])=TFR[4];(Subscript[g, 6])=TFR[4];(Subscript[h, 6])=TFR[4];(Subscript[e, 7])=TFR[4];(Subscript[f, 7])=TFR[4];(Subscript[g, 7])=TFR[4];(Subscript[h, 7])=TFR[4];(Subscript[e, 8])=TFR[4];(Subscript[f, 8])=TFR[4];(Subscript[g, 8])=TFR[4];(Subscript[h, 8])=TFR[4];

(Subscript[x, 1])=TFR[4] ;(Subscript[x, 2])=TFR[4] ;(Subscript[x, 3])=TFR[4] ;(Subscript[x, 4])=TFR[4];(Subscript[x, 5])=TFR[4] ;(Subscript[x, 6])=TFR[4] ;(Subscript[x, 7])=TFR[4];(Subscript[x, 8])=TFR[4] ;

(*  Let's look at at system with from 5 to 8 unknowns, omitting 1, 2 or 3 if there are just 7, 6 or 5 unknowns.
   *)

(Subscript[k, 1])=QuatDot[(Subscript[a, 1]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 1]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 1]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 1]),(Subscript[x, 4])]+QuatDot[(Subscript[e, 1]),(Subscript[x, 5])]+QuatDot[(Subscript[f, 1]),(Subscript[x, 6])]+QuatDot[(Subscript[g, 1]),(Subscript[x, 7])]+QuatDot[(Subscript[h, 1]),(Subscript[x, 8])];     (*  User must "ENTER" here.  *) 
(Subscript[k, 2])=QuatDot[(Subscript[a, 2]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 2]),(Subscript[x, 2])] +QuatDot[(Subscript[c, 2]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 2]),(Subscript[x, 4])]+QuatDot[(Subscript[e, 2]),(Subscript[x, 5])]+QuatDot[(Subscript[f, 2]),(Subscript[x, 6])] +QuatDot[(Subscript[g, 2]),(Subscript[x, 7])]+QuatDot[(Subscript[h, 2]),(Subscript[x, 8])];
(Subscript[k, 3])=QuatDot[(Subscript[a, 3]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 3]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 3]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 3]),(Subscript[x, 4])]+QuatDot[(Subscript[e, 3]),(Subscript[x, 5])]+QuatDot[(Subscript[f, 3]),(Subscript[x, 6])]+QuatDot[(Subscript[g, 3]),(Subscript[x, 7])]+QuatDot[(Subscript[h, 3]),(Subscript[x, 8])];
(Subscript[k, 4])=QuatDot[(Subscript[a, 4]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 4]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 4]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 4]),(Subscript[x, 4])]+QuatDot[(Subscript[e, 4]),(Subscript[x, 5])]+QuatDot[(Subscript[f, 4]),(Subscript[x, 6])]+QuatDot[(Subscript[g, 4]),(Subscript[x, 7])]+QuatDot[(Subscript[h, 4]),(Subscript[x, 8])];
(Subscript[k, 5])=QuatDot[(Subscript[a, 5]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 5]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 5]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 5]),(Subscript[x, 4])]+QuatDot[(Subscript[e, 5]),(Subscript[x, 5])]+QuatDot[(Subscript[f, 5]),(Subscript[x, 6])]+QuatDot[(Subscript[g, 5]),(Subscript[x, 7])]+QuatDot[(Subscript[h, 5]),(Subscript[x, 8])];
(Subscript[k, 6])=QuatDot[(Subscript[a, 6]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 6]),(Subscript[x, 2])] +QuatDot[(Subscript[c, 6]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 6]),(Subscript[x, 4])]+QuatDot[(Subscript[e, 6]),(Subscript[x, 5])]+QuatDot[(Subscript[f, 6]),(Subscript[x, 6])] +QuatDot[(Subscript[g, 6]),(Subscript[x, 7])]+QuatDot[(Subscript[h, 6]),(Subscript[x, 8])];
(Subscript[k, 7])=QuatDot[(Subscript[a, 7]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 7]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 7]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 7]),(Subscript[x, 4])]+QuatDot[(Subscript[e, 7]),(Subscript[x, 5])]+QuatDot[(Subscript[f, 7]),(Subscript[x, 6])]+QuatDot[(Subscript[g, 7]),(Subscript[x, 7])]+QuatDot[(Subscript[h, 7]),(Subscript[x, 8])];
(Subscript[k, 8])=QuatDot[(Subscript[a, 8]),(Subscript[x, 1])]+QuatDot[(Subscript[b, 8]),(Subscript[x, 2])]+QuatDot[(Subscript[c, 8]),(Subscript[x, 3])]+QuatDot[(Subscript[d, 8]),(Subscript[x, 4])]+QuatDot[(Subscript[e, 8]),(Subscript[x, 5])]+QuatDot[(Subscript[f, 8]),(Subscript[x, 6])]+QuatDot[(Subscript[g, 8]),(Subscript[x, 7])]+QuatDot[(Subscript[h, 8]),(Subscript[x, 8])];

A4=(Subscript[TransMatL, 4])[(Subscript[a, 1])];B4=(Subscript[TransMatL, 4])[(Subscript[b, 1])];C4=(Subscript[TransMatL, 4])[(Subscript[a, 2])];D4=(Subscript[TransMatL, 4])[(Subscript[b, 2])];A8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[c, 1])];B4=(Subscript[TransMatL, 4])[(Subscript[d, 1])];C4=(Subscript[TransMatL, 4])[(Subscript[c, 2])];D4=(Subscript[TransMatL, 4])[(Subscript[d, 2])];B8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[a, 3])];B4=(Subscript[TransMatL, 4])[(Subscript[b, 3])];C4=(Subscript[TransMatL, 4])[(Subscript[a, 4])];D4=(Subscript[TransMatL, 4])[(Subscript[b, 4])];C8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[c, 3])];B4=(Subscript[TransMatL, 4])[(Subscript[d, 3])];C4=(Subscript[TransMatL, 4])[(Subscript[c, 4])];D4=(Subscript[TransMatL, 4])[(Subscript[d, 4])];D8=BuildUp[4,A4,B4,C4,D4]

AA=BuildUp[8,A8,B8,C8,D8]

(*  At this point we have the first quarter built up for any PEP solution.  *)

A4=(Subscript[TransMatL, 4])[(Subscript[e, 1])];B4=(Subscript[TransMatL, 4])[(Subscript[f, 1])];C4=(Subscript[TransMatL, 4])[(Subscript[e, 2])];D4=(Subscript[TransMatL, 4])[(Subscript[f, 2])];A8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[g, 1])];B4=(Subscript[TransMatL, 4])[(Subscript[h, 1])];C4=(Subscript[TransMatL, 4])[(Subscript[g, 2])];D4=(Subscript[TransMatL, 4])[(Subscript[h, 2])];B8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[e, 3])];B4=(Subscript[TransMatL, 4])[(Subscript[f, 3])];C4=(Subscript[TransMatL, 4])[(Subscript[e, 4])];D4=(Subscript[TransMatL, 4])[(Subscript[f, 4])];C8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[g, 3])];B4=(Subscript[TransMatL, 4])[(Subscript[h, 3])];C4=(Subscript[TransMatL, 4])[(Subscript[g, 4])];D4=(Subscript[TransMatL, 4])[(Subscript[h, 4])];D8=BuildUp[4,A4,B4,C4,D4]

BB=BuildUp[8,A8,B8,C8,D8]

(*  We have the 2'nd quarter, BB.  *)

A4=(Subscript[TransMatL, 4])[(Subscript[a, 5])];B4=(Subscript[TransMatL, 4])[(Subscript[b, 5])];C4=(Subscript[TransMatL, 4])[(Subscript[a, 6])];D4=(Subscript[TransMatL, 4])[(Subscript[b, 6])];A8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[c, 5])];B4=(Subscript[TransMatL, 4])[(Subscript[d, 5])];C4=(Subscript[TransMatL, 4])[(Subscript[c, 6])];D4=(Subscript[TransMatL, 4])[(Subscript[d, 6])];B8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[a, 7])];B4=(Subscript[TransMatL, 4])[(Subscript[b, 7])];C4=(Subscript[TransMatL, 4])[(Subscript[a, 8])];D4=(Subscript[TransMatL, 4])[(Subscript[b, 8])];C8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[c, 7])];B4=(Subscript[TransMatL, 4])[(Subscript[d, 7])];C4=(Subscript[TransMatL, 4])[(Subscript[c, 8])];D4=(Subscript[TransMatL, 4])[(Subscript[d, 8])];D8=BuildUp[4,A4,B4,C4,D4]

CC=BuildUp[8,A8,B8,C8,D8]

(*  At this point we have the third quarter built up for any PEP solution.  *)

A4=(Subscript[TransMatL, 4])[(Subscript[e, 5])];B4=(Subscript[TransMatL, 4])[(Subscript[f, 5])];C4=(Subscript[TransMatL, 4])[(Subscript[e, 6])];D4=(Subscript[TransMatL, 4])[(Subscript[f, 6])];A8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[g, 5])];B4=(Subscript[TransMatL, 4])[(Subscript[h, 5])];C4=(Subscript[TransMatL, 4])[(Subscript[g, 6])];D4=(Subscript[TransMatL, 4])[(Subscript[h, 6])];B8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[e, 7])];B4=(Subscript[TransMatL, 4])[(Subscript[f, 7])];C4=(Subscript[TransMatL, 4])[(Subscript[e, 8])];D4=(Subscript[TransMatL, 4])[(Subscript[f, 8])];C8=BuildUp[4,A4,B4,C4,D4]

A4=(Subscript[TransMatL, 4])[(Subscript[g, 7])];B4=(Subscript[TransMatL, 4])[(Subscript[h, 7])];C4=(Subscript[TransMatL, 4])[(Subscript[g, 8])];D4=(Subscript[TransMatL, 4])[(Subscript[h, 8])];D8=BuildUp[4,A4,B4,C4,D4]

DD=BuildUp[8,A8,B8,C8,D8]

(*  We have the 4'th quarter, DD.  *)

st=AbsoluteTime[];(Subscript[PEP, 1])[AA,BB,CC,DD] ;AbsoluteTime[]-st  (*  This takes less than half a second.  *)

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[16]&&AA.X+BB.Z==ZeroMat[16]&&CC.W+DD.Y==ZeroMat[16]&&CC.X+DD.Z==IdentityMatrix[16]

W.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+X.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]==Join[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4])]

Y.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+Z.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]

Join[(Subscript[x, 5]),(Subscript[x, 6]),(Subscript[x, 7]),(Subscript[x, 8])]

{Dimensions[W ],((Length[W])^2)}    (*  The matrix inverse upper left quarter, W, takes up 3+ screens worth, if you're curious.  *)

(*  ----------------------------------  fin: 8x8's Quaternion Systems ------------------------------------------------  *)

(*  To illustrate the solution of a 4x4 system, we randomly select 16 octonion constants 
  with +/- digit components.  *)

SeedRandom[9797];(Subscript[a, 1])=TFR[8];(Subscript[a, 2])=TFR[8];(Subscript[b, 1])=TFR[8];(Subscript[b, 2])=TFR[8];(Subscript[c, 1])=TFR[8];(Subscript[c, 2])=TFR[8];(Subscript[d, 1])=TFR[8];(Subscript[d, 2])=TFR[8];(Subscript[a, 3])=TFR[8];(Subscript[a, 4])=TFR[8];(Subscript[b, 3])=TFR[8];(Subscript[b, 4])=TFR[8];(Subscript[c, 3])=TFR[8];(Subscript[c, 4])=TFR[8];(Subscript[d, 3])=TFR[8];(Subscript[d, 4])=TFR[8];
(*  User must "ENTER" here. You may (of course) change the seed number above to generate new systems.  *)

(Subscript[a, 1])

(*   As is often done in textbook examples, randomly select the answers so that the solution will avoid 
those ugly fractions.  *)

(*  "Answers" for (Subscript[x, 1]), ..., (Subscript[x, 4])  *)(*  User must "ENTER" below.  *)

(Subscript[x, 1])=TFR[8] ;(Subscript[x, 2])=TFR[8] ;(Subscript[x, 3])=TFR[8] ;(Subscript[x, 4])=TFR[8];

(*  Compute the rhs constants for multiplication from the left.  *)

(Subscript[k, 1])=OctDot[(Subscript[a, 1]),(Subscript[x, 1])]+OctDot[(Subscript[b, 1]),(Subscript[x, 2])]+OctDot[(Subscript[c, 1]),(Subscript[x, 3])]+OctDot[(Subscript[d, 1]),(Subscript[x, 4])];(Subscript[k, 2])=OctDot[(Subscript[a, 2]),(Subscript[x, 1])]+OctDot[(Subscript[b, 2]),(Subscript[x, 2])] +OctDot[(Subscript[c, 2]),(Subscript[x, 3])]+OctDot[(Subscript[d, 2]),(Subscript[x, 4])];(Subscript[k, 3])=OctDot[(Subscript[a, 3]),(Subscript[x, 1])]+OctDot[(Subscript[b, 3]),(Subscript[x, 2])]+OctDot[(Subscript[c, 3]),(Subscript[x, 3])]+OctDot[(Subscript[d, 3]),(Subscript[x, 4])];(Subscript[k, 4])=OctDot[(Subscript[a, 4]),(Subscript[x, 1])]+OctDot[(Subscript[b, 4]),(Subscript[x, 2])]+OctDot[(Subscript[c, 4]),(Subscript[x, 3])]+OctDot[(Subscript[d, 4]),(Subscript[x, 4])]             (*  User must "ENTER" here.  *)

(*  We illustrate solving the octonion system:
	 (Subscript[a, 1])(Subscript[x, 1]) + (Subscript[b, 1])(Subscript[x, 2]) + (Subscript[c, 1])(Subscript[x, 3]) + (Subscript[d, 1])(Subscript[x, 4]) = (Subscript[k, 1]),
		(Subscript[a, 2])(Subscript[x, 1]) + (Subscript[b, 2])(Subscript[x, 2]) + (Subscript[c, 2])(Subscript[x, 3]) + (Subscript[d, 2])(Subscript[x, 4]) = (Subscript[k, 2]),
		(Subscript[a, 3])(Subscript[x, 1]) + (Subscript[b, 3])(Subscript[x, 2]) + (Subscript[c, 3])(Subscript[x, 3]) + (Subscript[d, 3])(Subscript[x, 4]) = (Subscript[k, 3]),
		(Subscript[a, 4])(Subscript[x, 1]) + (Subscript[b, 4])(Subscript[x, 2]) + (Subscript[c, 4])(Subscript[x, 3]) + (Subscript[d, 4])(Subscript[x, 4]) = (Subscript[k, 4]),		
  *)

A8=(Subscript[TransMatL, 8])[(Subscript[a, 1])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 1])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 2])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 2])];AA=BuildUp[8,A8,B8,C8,D8]

A8=(Subscript[TransMatL, 8])[(Subscript[c, 1])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 1])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 2])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 2])];BB=BuildUp[8,A8,B8,C8,D8]

A8=(Subscript[TransMatL, 8])[(Subscript[a, 3])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 3])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 4])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 4])];CC=BuildUp[8,A8,B8,C8,D8]

A8=(Subscript[TransMatL, 8])[(Subscript[c, 3])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 3])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 4])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 4])];DD=BuildUp[8,A8,B8,C8,D8]

(Subscript[PEP, 2])[AA,BB,CC,DD]

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[16]&&AA.X+BB.Z==ZeroMat[16]&&CC.W+DD.Y==ZeroMat[16]&&CC.X+DD.Z==IdentityMatrix[16]

(Subscript[x, 1])

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]

(Subscript[x, 2])

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 1]),(Subscript[x, 2])]&&Y.Join[(Subscript[k, 1]),(Subscript[k, 2])]+Z.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 3]),(Subscript[x, 4])]

(*  Compute the rhs constants for multiplication from the right.  *)

(Subscript[k, 1])=OctDot[(Subscript[x, 1]),(Subscript[a, 1])]+OctDot[(Subscript[x, 2]),(Subscript[b, 1])]+OctDot[(Subscript[x, 3]),(Subscript[c, 1])]+OctDot[(Subscript[x, 4]),(Subscript[d, 1])];(Subscript[k, 2])=OctDot[(Subscript[x, 1]),(Subscript[a, 2])]+OctDot[(Subscript[x, 2]),(Subscript[b, 2])] +OctDot[(Subscript[x, 3]),(Subscript[c, 2])]+OctDot[(Subscript[x, 4]),(Subscript[d, 2])];(Subscript[k, 3])=OctDot[(Subscript[x, 1]),(Subscript[a, 3])]+OctDot[(Subscript[x, 2]),(Subscript[b, 3])]+OctDot[(Subscript[x, 3]),(Subscript[c, 3])]+OctDot[(Subscript[x, 4]),(Subscript[d, 3])];(Subscript[k, 4])=OctDot[(Subscript[x, 1]),(Subscript[a, 4])]+OctDot[(Subscript[x, 2]),(Subscript[b, 4])]+OctDot[(Subscript[x, 3]),(Subscript[c, 4])]+OctDot[(Subscript[x, 4]),(Subscript[d, 4])]             (*  User must "ENTER" here.  *)

(*  We illustrate solving the octonion system:
	 (Subscript[x, 1])(Subscript[a, 1]) + (Subscript[x, 2])(Subscript[b, 1]) + (Subscript[x, 3])(Subscript[c, 1]) + (Subscript[x, 4])(Subscript[d, 1]) = (Subscript[k, 1]),
		(Subscript[x, 1])(Subscript[a, 2]) + (Subscript[x, 2])(Subscript[b, 2]) + (Subscript[x, 3])(Subscript[c, 2]) + (Subscript[x, 4])(Subscript[d, 2]) = (Subscript[k, 2]),
		(Subscript[x, 1])(Subscript[a, 3]) + (Subscript[x, 2])(Subscript[b, 3]) + (Subscript[x, 3])(Subscript[c, 3]) + (Subscript[x, 4])(Subscript[d, 3]) = (Subscript[k, 3]),
		(Subscript[x, 1])(Subscript[a, 4]) + (Subscript[x, 2])(Subscript[b, 4]) + (Subscript[x, 3])(Subscript[c, 4]) + (Subscript[x, 4])(Subscript[d, 4]) = (Subscript[k, 4]),		
  *)

A8=(Subscript[TransMatR, 8])[(Subscript[a, 1])];B8=(Subscript[TransMatR, 8])[(Subscript[b, 1])];C8=(Subscript[TransMatR, 8])[(Subscript[a, 2])];D8=(Subscript[TransMatR, 8])[(Subscript[b, 2])];AA=BuildUp[8,A8,B8,C8,D8]

A8=(Subscript[TransMatR, 8])[(Subscript[c, 1])];B8=(Subscript[TransMatR, 8])[(Subscript[d, 1])];C8=(Subscript[TransMatR, 8])[(Subscript[c, 2])];D8=(Subscript[TransMatR, 8])[(Subscript[d, 2])];BB=BuildUp[8,A8,B8,C8,D8]

A8=(Subscript[TransMatR, 8])[(Subscript[a, 3])];B8=(Subscript[TransMatR, 8])[(Subscript[b, 3])];C8=(Subscript[TransMatR, 8])[(Subscript[a, 4])];D8=(Subscript[TransMatR, 8])[(Subscript[b, 4])];CC=BuildUp[8,A8,B8,C8,D8]

A8=(Subscript[TransMatR, 8])[(Subscript[c, 3])];B8=(Subscript[TransMatR, 8])[(Subscript[d, 3])];C8=(Subscript[TransMatR, 8])[(Subscript[c, 4])];D8=(Subscript[TransMatR, 8])[(Subscript[d, 4])];DD=BuildUp[8,A8,B8,C8,D8]

(Subscript[PEP, 3])[AA,BB,CC,DD]

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[16]&&AA.X+BB.Z==ZeroMat[16]&&CC.W+DD.Y==ZeroMat[16]&&CC.X+DD.Z==IdentityMatrix[16]

(Subscript[x, 1])

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]

(Subscript[x, 2])

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 1]),(Subscript[x, 2])]&&Y.Join[(Subscript[k, 1]),(Subscript[k, 2])]+Z.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 3]),(Subscript[x, 4])]

(*  We select 64 constants as coefficients of an 8x8 system of octonions.  *)

SeedRandom[9797];(Subscript[a, 1])=TFR[8];(Subscript[a, 2])=TFR[8];(Subscript[b, 1])=TFR[8];(Subscript[b, 2])=TFR[8];(Subscript[c, 1])=TFR[8];(Subscript[c, 2])=TFR[8];(Subscript[d, 1])=TFR[8];(Subscript[d, 2])=TFR[8];(Subscript[a, 3])=TFR[8];(Subscript[a, 4])=TFR[8];(Subscript[b, 3])=TFR[8];(Subscript[b, 4])=TFR[8];(Subscript[c, 3])=TFR[8];(Subscript[c, 4])=TFR[8];(Subscript[d, 3])=TFR[8];(Subscript[d, 4])=TFR[8];(Subscript[e, 1])=TFR[8];(Subscript[f, 1])=TFR[8];(Subscript[g, 1])=TFR[8];(Subscript[h, 1])=TFR[8];(Subscript[e, 2])=TFR[8];(Subscript[f, 2])=TFR[8];(Subscript[g, 2])=TFR[8];(Subscript[h, 2])=TFR[8];(Subscript[e, 3])=TFR[8];(Subscript[f, 3])=TFR[8];(Subscript[g, 3])=TFR[8];(Subscript[h, 3])=TFR[8];(Subscript[e, 4])=TFR[8];(Subscript[f, 4])=TFR[8];(Subscript[g, 4])=TFR[8];(Subscript[h, 4])=TFR[8];(Subscript[e, 5])=TFR[8];(Subscript[f, 5])=TFR[8];(Subscript[g, 5])=TFR[8];(Subscript[h, 5])=TFR[8];(Subscript[e, 6])=TFR[8];(Subscript[f, 6])=TFR[8];(Subscript[g, 6])=TFR[8];(Subscript[h, 6])=TFR[8];(Subscript[e, 7])=TFR[8];(Subscript[f, 7])=TFR[8];(Subscript[g, 7])=TFR[8];(Subscript[h, 7])=TFR[8];(Subscript[e, 8])=TFR[8];(Subscript[f, 8])=TFR[8];(Subscript[g, 8])=TFR[8];(Subscript[h, 8])=TFR[8];(Subscript[a, 5])=TFR[8];(Subscript[b, 5])=TFR[8];(Subscript[c, 5])=TFR[8];(Subscript[d, 5])=TFR[8];(Subscript[a, 6])=TFR[8];(Subscript[b, 6])=TFR[8];(Subscript[c, 6])=TFR[8];(Subscript[d, 6])=TFR[8];(Subscript[a, 7])=TFR[8];(Subscript[b, 7])=TFR[8];(Subscript[c, 7])=TFR[8];(Subscript[d, 7])=TFR[8];(Subscript[a, 8])=TFR[8];(Subscript[b, 8])=TFR[8];(Subscript[c, 8])=TFR[8];(Subscript[d, 8])=TFR[8];
(*  User must "ENTER" here. You may (of course) change the seed number above to generate new systems.  *)

(Subscript[x, 8])=TFR[8];(Subscript[x, 7])=TFR[8];(Subscript[x, 6])=TFR[8];(Subscript[x, 5])=TFR[8];(Subscript[x, 4])=TFR[8];(Subscript[x, 3])=TFR[8];(Subscript[x, 2])=TFR[8];(Subscript[x, 1])=TFR[8];

(*  We have the system with constant multiplying from the right side.  *)

(*  User must "ENTER" here.  *) (Subscript[k, 1])=OctDot[(Subscript[x, 1]),(Subscript[a, 1])]+OctDot[(Subscript[x, 2]),(Subscript[b, 1])]+OctDot[(Subscript[x, 3]),(Subscript[c, 1])]+OctDot[(Subscript[x, 4]),(Subscript[d, 1])]+OctDot[(Subscript[x, 5]),(Subscript[e, 1])]+OctDot[(Subscript[x, 6]),(Subscript[f, 1])]+OctDot[(Subscript[x, 7]),(Subscript[g, 1])]+OctDot[(Subscript[x, 8]),(Subscript[h, 1])];     
(Subscript[k, 2])=OctDot[(Subscript[x, 1]),(Subscript[a, 2])]+OctDot[(Subscript[x, 2]),(Subscript[b, 2])] +OctDot[(Subscript[x, 3]),(Subscript[c, 2])]+OctDot[(Subscript[x, 4]),(Subscript[d, 2])]+OctDot[(Subscript[x, 5]),(Subscript[e, 2])]+OctDot[(Subscript[x, 6]),(Subscript[f, 2])]+OctDot[(Subscript[x, 7]),(Subscript[g, 2])]+OctDot[(Subscript[x, 8]),(Subscript[h, 2])];
(Subscript[k, 3])=OctDot[(Subscript[x, 1]),(Subscript[a, 3])]+OctDot[(Subscript[x, 2]),(Subscript[b, 3])]+OctDot[(Subscript[x, 3]),(Subscript[c, 3])]+OctDot[(Subscript[x, 4]),(Subscript[d, 3])]+OctDot[(Subscript[x, 5]),(Subscript[e, 3])]+OctDot[(Subscript[x, 6]),(Subscript[f, 3])]+OctDot[(Subscript[x, 7]),(Subscript[g, 3])]+OctDot[(Subscript[x, 8]),(Subscript[h, 3])];
(Subscript[k, 4])=OctDot[(Subscript[x, 1]),(Subscript[a, 4])]+OctDot[(Subscript[x, 2]),(Subscript[b, 4])]+OctDot[(Subscript[x, 3]),(Subscript[c, 4])]+OctDot[(Subscript[x, 4]),(Subscript[d, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 4])]+OctDot[(Subscript[x, 6]),(Subscript[f, 4])]+OctDot[(Subscript[x, 7]),(Subscript[g, 4])]+OctDot[(Subscript[x, 8]),(Subscript[h, 4])];
(Subscript[k, 5])=OctDot[(Subscript[x, 1]),(Subscript[a, 5])]+OctDot[(Subscript[x, 2]),(Subscript[b, 5])]+OctDot[(Subscript[x, 3]),(Subscript[c, 5])]+OctDot[(Subscript[x, 4]),(Subscript[d, 5])]+OctDot[(Subscript[x, 5]),(Subscript[e, 5])]+OctDot[(Subscript[x, 6]),(Subscript[f, 5])]+OctDot[(Subscript[x, 7]),(Subscript[g, 5])]+OctDot[(Subscript[x, 8]),(Subscript[h, 5])];
(Subscript[k, 6])=OctDot[(Subscript[x, 1]),(Subscript[a, 6])]+OctDot[(Subscript[x, 2]),(Subscript[b, 6])] +OctDot[(Subscript[x, 3]),(Subscript[c, 6])]+OctDot[(Subscript[x, 4]),(Subscript[d, 6])]+OctDot[(Subscript[x, 5]),(Subscript[e, 6])]+OctDot[(Subscript[x, 6]),(Subscript[f, 6])] +OctDot[(Subscript[x, 7]),(Subscript[g, 6])]+OctDot[(Subscript[x, 8]),(Subscript[h, 6])];
(Subscript[k, 7])=OctDot[(Subscript[x, 1]),(Subscript[a, 7])]+OctDot[(Subscript[x, 2]),(Subscript[b, 7])]+OctDot[(Subscript[x, 3]),(Subscript[c, 7])]+OctDot[(Subscript[x, 4]),(Subscript[d, 7])]+OctDot[(Subscript[x, 5]),(Subscript[e, 7])]+OctDot[(Subscript[x, 6]),(Subscript[f, 7])]+OctDot[(Subscript[x, 7]),(Subscript[g, 7])]+OctDot[(Subscript[x, 8]),(Subscript[h, 7])];
(Subscript[k, 8])=OctDot[(Subscript[x, 1]),(Subscript[a, 8])]+OctDot[(Subscript[x, 2]),(Subscript[b, 8])]+OctDot[(Subscript[x, 3]),(Subscript[c, 8])]+OctDot[(Subscript[x, 4]),(Subscript[d, 8])]+OctDot[(Subscript[x, 5]),(Subscript[e, 8])]+OctDot[(Subscript[x, 6]),(Subscript[f, 8])]+OctDot[(Subscript[x, 7]),(Subscript[g, 8])]+OctDot[(Subscript[x, 8]),(Subscript[h, 8])];

{(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4]),(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])}

A8=(Subscript[TransMatR, 8])[(Subscript[a, 1])];B8=(Subscript[TransMatR, 8])[(Subscript[b, 1])];C8=(Subscript[TransMatR, 8])[(Subscript[a, 2])];D8=(Subscript[TransMatR, 8])[(Subscript[b, 2])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[c, 1])];B8=(Subscript[TransMatR, 8])[(Subscript[d, 1])];C8=(Subscript[TransMatR, 8])[(Subscript[c, 2])];D8=(Subscript[TransMatR, 8])[(Subscript[d, 2])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[a, 3])];B8=(Subscript[TransMatR, 8])[(Subscript[b, 3])];C8=(Subscript[TransMatR, 8])[(Subscript[a, 4])];D8=(Subscript[TransMatR, 8])[(Subscript[b, 4])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[c, 3])];B8=(Subscript[TransMatR, 8])[(Subscript[d, 3])];C8=(Subscript[TransMatR, 8])[(Subscript[c, 4])];D8=(Subscript[TransMatR, 8])[(Subscript[d, 4])];D16=BuildUp[8,A8,B8,C8,D8];AA=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatR, 8])[(Subscript[e, 1])];B8=(Subscript[TransMatR, 8])[(Subscript[f, 1])];C8=(Subscript[TransMatR, 8])[(Subscript[e, 2])];D8=(Subscript[TransMatR, 8])[(Subscript[f, 2])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[g, 1])];B8=(Subscript[TransMatR, 8])[(Subscript[h, 1])];C8=(Subscript[TransMatR, 8])[(Subscript[g, 2])];D8=(Subscript[TransMatR, 8])[(Subscript[h, 2])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[e, 3])];B8=(Subscript[TransMatR, 8])[(Subscript[f, 3])];C8=(Subscript[TransMatR, 8])[(Subscript[e, 4])];D8=(Subscript[TransMatR, 8])[(Subscript[f, 4])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[g, 3])];B8=(Subscript[TransMatR, 8])[(Subscript[h, 3])];C8=(Subscript[TransMatR, 8])[(Subscript[g, 4])];D8=(Subscript[TransMatR, 8])[(Subscript[h, 4])];D16=BuildUp[8,A8,B8,C8,D8];BB=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatR, 8])[(Subscript[a, 5])];B8=(Subscript[TransMatR, 8])[(Subscript[b, 5])];C8=(Subscript[TransMatR, 8])[(Subscript[a, 6])];D8=(Subscript[TransMatR, 8])[(Subscript[b, 6])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[c, 5])];B8=(Subscript[TransMatR, 8])[(Subscript[d, 5])];C8=(Subscript[TransMatR, 8])[(Subscript[c, 6])];D8=(Subscript[TransMatR, 8])[(Subscript[d, 6])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[a, 7])];B8=(Subscript[TransMatR, 8])[(Subscript[b, 7])];C8=(Subscript[TransMatR, 8])[(Subscript[a, 8])];D8=(Subscript[TransMatR, 8])[(Subscript[b, 8])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[c, 7])];B8=(Subscript[TransMatR, 8])[(Subscript[d, 7])];C8=(Subscript[TransMatR, 8])[(Subscript[c, 8])];D8=(Subscript[TransMatR, 8])[(Subscript[d, 8])];D16=BuildUp[8,A8,B8,C8,D8];CC=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatR, 8])[(Subscript[e, 5])];B8=(Subscript[TransMatR, 8])[(Subscript[f, 5])];C8=(Subscript[TransMatR, 8])[(Subscript[e, 6])];D8=(Subscript[TransMatR, 8])[(Subscript[f, 6])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[g, 5])];B8=(Subscript[TransMatR, 8])[(Subscript[h, 5])];C8=(Subscript[TransMatR, 8])[(Subscript[g, 6])];D8=(Subscript[TransMatR, 8])[(Subscript[h, 6])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[e, 7])];B8=(Subscript[TransMatR, 8])[(Subscript[f, 7])];C8=(Subscript[TransMatR, 8])[(Subscript[e, 8])];D8=(Subscript[TransMatR, 8])[(Subscript[f, 8])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[g, 7])];B8=(Subscript[TransMatR, 8])[(Subscript[h, 7])];C8=(Subscript[TransMatR, 8])[(Subscript[g, 8])];D8=(Subscript[TransMatR, 8])[(Subscript[h, 8])];D16=BuildUp[8,A8,B8,C8,D8];DD=BuildUp[16,A16,B16,C16,D16];

st=AbsoluteTime[];(Subscript[PEP, 2])[AA,BB,CC,DD] ;AbsoluteTime[]-st  (*  This takes less than half a dozen seconds.  *)

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[32]&&AA.X+BB.Z==ZeroMat[32]&&CC.W+DD.Y==ZeroMat[32]&&CC.X+DD.Z==IdentityMatrix[32]

{(Subscript[x, 1]),(Subscript[x, 2])}

Floor[N[W.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+X.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]]]

{(Subscript[x, 3]),(Subscript[x, 4])}

Y.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+Z.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]==Join[(Subscript[x, 5]),(Subscript[x, 6]),(Subscript[x, 7]),(Subscript[x, 8])]

(*  We can illustrate multultiplication from the left for the abcd coefficients, leaving multiplication from the right 
  for the efgh coefficients when computing the constants.  *)

(*  User must "ENTER" here.  *) (Subscript[k, 1])=OctDot[(Subscript[a, 1]),(Subscript[x, 1])]+OctDot[(Subscript[b, 1]),(Subscript[x, 2])]+OctDot[(Subscript[c, 1]),(Subscript[x, 3])]+OctDot[(Subscript[d, 1]),(Subscript[x, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 1])]+OctDot[(Subscript[x, 6]),(Subscript[f, 1])]+OctDot[(Subscript[x, 7]),(Subscript[g, 1])]+OctDot[(Subscript[x, 8]),(Subscript[h, 1])];     
(Subscript[k, 2])=OctDot[(Subscript[a, 2]),(Subscript[x, 1])]+OctDot[(Subscript[b, 2]),(Subscript[x, 2])]+OctDot[(Subscript[c, 2]),(Subscript[x, 3])]+OctDot[(Subscript[d, 2]),(Subscript[x, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 2])]+OctDot[(Subscript[x, 6]),(Subscript[f, 2])]+OctDot[(Subscript[x, 7]),(Subscript[g, 2])]+OctDot[(Subscript[x, 8]),(Subscript[h, 2])];
(Subscript[k, 3])=OctDot[(Subscript[a, 3]),(Subscript[x, 1])]+OctDot[(Subscript[b, 3]),(Subscript[x, 2])]+OctDot[(Subscript[c, 3]),(Subscript[x, 3])]+OctDot[(Subscript[d, 3]),(Subscript[x, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 3])]+OctDot[(Subscript[x, 6]),(Subscript[f, 3])]+OctDot[(Subscript[x, 7]),(Subscript[g, 3])]+OctDot[(Subscript[x, 8]),(Subscript[h, 3])];
(Subscript[k, 4])=OctDot[(Subscript[a, 4]),(Subscript[x, 1])]+OctDot[(Subscript[b, 4]),(Subscript[x, 2])]+OctDot[(Subscript[c, 4]),(Subscript[x, 3])]+OctDot[(Subscript[d, 4]),(Subscript[x, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 4])]+OctDot[(Subscript[x, 6]),(Subscript[f, 4])]+OctDot[(Subscript[x, 7]),(Subscript[g, 4])]+OctDot[(Subscript[x, 8]),(Subscript[h, 4])];
(Subscript[k, 5])=OctDot[(Subscript[a, 5]),(Subscript[x, 1])]+OctDot[(Subscript[b, 5]),(Subscript[x, 2])]+OctDot[(Subscript[c, 5]),(Subscript[x, 3])]+OctDot[(Subscript[d, 5]),(Subscript[x, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 5])]+OctDot[(Subscript[x, 6]),(Subscript[f, 5])]+OctDot[(Subscript[x, 7]),(Subscript[g, 5])]+OctDot[(Subscript[x, 8]),(Subscript[h, 5])];
(Subscript[k, 6])=OctDot[(Subscript[a, 6]),(Subscript[x, 1])]+OctDot[(Subscript[b, 6]),(Subscript[x, 2])]+OctDot[(Subscript[c, 6]),(Subscript[x, 3])]+OctDot[(Subscript[d, 6]),(Subscript[x, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 6])]+OctDot[(Subscript[x, 6]),(Subscript[f, 6])]+OctDot[(Subscript[x, 7]),(Subscript[g, 6])]+OctDot[(Subscript[x, 8]),(Subscript[h, 6])];
(Subscript[k, 7])=OctDot[(Subscript[a, 7]),(Subscript[x, 1])]+OctDot[(Subscript[b, 7]),(Subscript[x, 2])]+OctDot[(Subscript[c, 7]),(Subscript[x, 3])]+OctDot[(Subscript[d, 7]),(Subscript[x, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 7])]+OctDot[(Subscript[x, 6]),(Subscript[f, 7])]+OctDot[(Subscript[x, 7]),(Subscript[g, 7])]+OctDot[(Subscript[x, 8]),(Subscript[h, 7])];
(Subscript[k, 8])=OctDot[(Subscript[a, 8]),(Subscript[x, 1])]+OctDot[(Subscript[b, 8]),(Subscript[x, 2])]+OctDot[(Subscript[c, 8]),(Subscript[x, 3])]+OctDot[(Subscript[d, 8]),(Subscript[x, 4])]+OctDot[(Subscript[x, 5]),(Subscript[e, 8])]+OctDot[(Subscript[x, 6]),(Subscript[f, 8])]+OctDot[(Subscript[x, 7]),(Subscript[g, 8])]+OctDot[(Subscript[x, 8]),(Subscript[h, 8])];

A8=(Subscript[TransMatL, 8])[(Subscript[a, 1])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 1])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 2])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 2])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[c, 1])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 1])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 2])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 2])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[a, 3])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 3])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 4])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 4])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[c, 3])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 3])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 4])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 4])];D16=BuildUp[8,A8,B8,C8,D8];AA=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatR, 8])[(Subscript[e, 1])];B8=(Subscript[TransMatR, 8])[(Subscript[f, 1])];C8=(Subscript[TransMatR, 8])[(Subscript[e, 2])];D8=(Subscript[TransMatR, 8])[(Subscript[f, 2])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[g, 1])];B8=(Subscript[TransMatR, 8])[(Subscript[h, 1])];C8=(Subscript[TransMatR, 8])[(Subscript[g, 2])];D8=(Subscript[TransMatR, 8])[(Subscript[h, 2])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[e, 3])];B8=(Subscript[TransMatR, 8])[(Subscript[f, 3])];C8=(Subscript[TransMatR, 8])[(Subscript[e, 4])];D8=(Subscript[TransMatR, 8])[(Subscript[f, 4])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[g, 3])];B8=(Subscript[TransMatR, 8])[(Subscript[h, 3])];C8=(Subscript[TransMatR, 8])[(Subscript[g, 4])];D8=(Subscript[TransMatR, 8])[(Subscript[h, 4])];D16=BuildUp[8,A8,B8,C8,D8];BB=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatL, 8])[(Subscript[a, 5])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 5])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 6])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 6])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[c, 5])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 5])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 6])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 6])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[a, 7])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 7])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 8])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 8])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[c, 7])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 7])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 8])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 8])];D16=BuildUp[8,A8,B8,C8,D8];CC=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatR, 8])[(Subscript[e, 5])];B8=(Subscript[TransMatR, 8])[(Subscript[f, 5])];C8=(Subscript[TransMatR, 8])[(Subscript[e, 6])];D8=(Subscript[TransMatR, 8])[(Subscript[f, 6])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[g, 5])];B8=(Subscript[TransMatR, 8])[(Subscript[h, 5])];C8=(Subscript[TransMatR, 8])[(Subscript[g, 6])];D8=(Subscript[TransMatR, 8])[(Subscript[h, 6])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[e, 7])];B8=(Subscript[TransMatR, 8])[(Subscript[f, 7])];C8=(Subscript[TransMatR, 8])[(Subscript[e, 8])];D8=(Subscript[TransMatR, 8])[(Subscript[f, 8])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatR, 8])[(Subscript[g, 7])];B8=(Subscript[TransMatR, 8])[(Subscript[h, 7])];C8=(Subscript[TransMatR, 8])[(Subscript[g, 8])];D8=(Subscript[TransMatR, 8])[(Subscript[h, 8])];D16=BuildUp[8,A8,B8,C8,D8];DD=BuildUp[16,A16,B16,C16,D16];

st=AbsoluteTime[];(Subscript[PEP, 3])[AA,BB,CC,DD] ;AbsoluteTime[]-st  (*  This takes less than half a dozen seconds.  *)

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[32]&&AA.X+BB.Z==ZeroMat[32]&&CC.W+DD.Y==ZeroMat[32]&&CC.X+DD.Z==IdentityMatrix[32]

{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4])}

Floor[N[W.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+X.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]]]

{-1,-9,4,8,-8,0,-8,6,      -9,2,-9,-3,0,-2,-4,-5,      -3,-7,1,-8,6,-9,5,2,    -8,-3,2,-3,-3,1,2,-7}

Y.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+Z.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]==Join[(Subscript[x, 5]),(Subscript[x, 6]),(Subscript[x, 7]),(Subscript[x, 8])]

(*  Now let all 64 coefficients be multiplying from the left  *)

(*  User must "ENTER" here.  *) (Subscript[k, 1])=OctDot[(Subscript[a, 1]),(Subscript[x, 1])]+OctDot[(Subscript[b, 1]),(Subscript[x, 2])]+OctDot[(Subscript[c, 1]),(Subscript[x, 3])]+OctDot[(Subscript[d, 1]),(Subscript[x, 4])]+OctDot[(Subscript[e, 1]),(Subscript[x, 5])]+OctDot[(Subscript[f, 1]),(Subscript[x, 6])]+OctDot[(Subscript[g, 1]),(Subscript[x, 7])]+OctDot[(Subscript[h, 1]),(Subscript[x, 8])];     
(Subscript[k, 2])=OctDot[(Subscript[a, 2]),(Subscript[x, 1])]+OctDot[(Subscript[b, 2]),(Subscript[x, 2])]+OctDot[(Subscript[c, 2]),(Subscript[x, 3])]+OctDot[(Subscript[d, 2]),(Subscript[x, 4])]+OctDot[(Subscript[e, 2]),(Subscript[x, 5])]+OctDot[(Subscript[f, 2]),(Subscript[x, 6])]+OctDot[(Subscript[g, 2]),(Subscript[x, 7])]+OctDot[(Subscript[h, 2]),(Subscript[x, 8])];
(Subscript[k, 3])=OctDot[(Subscript[a, 3]),(Subscript[x, 1])]+OctDot[(Subscript[b, 3]),(Subscript[x, 2])]+OctDot[(Subscript[c, 3]),(Subscript[x, 3])]+OctDot[(Subscript[d, 3]),(Subscript[x, 4])]+OctDot[(Subscript[e, 3]),(Subscript[x, 5])]+OctDot[(Subscript[f, 3]),(Subscript[x, 6])]+OctDot[(Subscript[g, 3]),(Subscript[x, 7])]+OctDot[(Subscript[h, 3]),(Subscript[x, 8])];
(Subscript[k, 4])=OctDot[(Subscript[a, 4]),(Subscript[x, 1])]+OctDot[(Subscript[b, 4]),(Subscript[x, 2])]+OctDot[(Subscript[c, 4]),(Subscript[x, 3])]+OctDot[(Subscript[d, 4]),(Subscript[x, 4])]+OctDot[(Subscript[e, 4]),(Subscript[x, 5])]+OctDot[(Subscript[f, 4]),(Subscript[x, 6])]+OctDot[(Subscript[g, 4]),(Subscript[x, 7])]+OctDot[(Subscript[h, 4]),(Subscript[x, 8])];
(Subscript[k, 5])=OctDot[(Subscript[a, 5]),(Subscript[x, 1])]+OctDot[(Subscript[b, 5]),(Subscript[x, 2])]+OctDot[(Subscript[c, 5]),(Subscript[x, 3])]+OctDot[(Subscript[d, 5]),(Subscript[x, 4])]+OctDot[(Subscript[e, 5]),(Subscript[x, 5])]+OctDot[(Subscript[f, 5]),(Subscript[x, 6])]+OctDot[(Subscript[g, 5]),(Subscript[x, 7])]+OctDot[(Subscript[h, 5]),(Subscript[x, 8])];
(Subscript[k, 6])=OctDot[(Subscript[a, 6]),(Subscript[x, 1])]+OctDot[(Subscript[b, 6]),(Subscript[x, 2])]+OctDot[(Subscript[c, 6]),(Subscript[x, 3])]+OctDot[(Subscript[d, 6]),(Subscript[x, 4])]+OctDot[(Subscript[e, 6]),(Subscript[x, 5])]+OctDot[(Subscript[f, 6]),(Subscript[x, 6])]+OctDot[(Subscript[g, 6]),(Subscript[x, 7])]+OctDot[(Subscript[h, 6]),(Subscript[x, 8])];
(Subscript[k, 7])=OctDot[(Subscript[a, 7]),(Subscript[x, 1])]+OctDot[(Subscript[b, 7]),(Subscript[x, 2])]+OctDot[(Subscript[c, 7]),(Subscript[x, 3])]+OctDot[(Subscript[d, 7]),(Subscript[x, 4])]+OctDot[(Subscript[e, 7]),(Subscript[x, 5])]+OctDot[(Subscript[f, 7]),(Subscript[x, 6])]+OctDot[(Subscript[g, 7]),(Subscript[x, 7])]+OctDot[(Subscript[h, 7]),(Subscript[x, 8])];
(Subscript[k, 8])=OctDot[(Subscript[a, 8]),(Subscript[x, 1])]+OctDot[(Subscript[b, 8]),(Subscript[x, 2])]+OctDot[(Subscript[c, 8]),(Subscript[x, 3])]+OctDot[(Subscript[d, 8]),(Subscript[x, 4])]+OctDot[(Subscript[e, 8]),(Subscript[x, 5])]+OctDot[(Subscript[f, 8]),(Subscript[x, 6])]+OctDot[(Subscript[g, 8]),(Subscript[x, 7])]+OctDot[(Subscript[h, 8]),(Subscript[x, 8])];

A8=(Subscript[TransMatL, 8])[(Subscript[a, 1])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 1])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 2])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 2])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[c, 1])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 1])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 2])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 2])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[a, 3])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 3])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 4])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 4])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[c, 3])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 3])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 4])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 4])];D16=BuildUp[8,A8,B8,C8,D8];AA=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatL, 8])[(Subscript[e, 1])];B8=(Subscript[TransMatL, 8])[(Subscript[f, 1])];C8=(Subscript[TransMatL, 8])[(Subscript[e, 2])];D8=(Subscript[TransMatL, 8])[(Subscript[f, 2])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[g, 1])];B8=(Subscript[TransMatL, 8])[(Subscript[h, 1])];C8=(Subscript[TransMatL, 8])[(Subscript[g, 2])];D8=(Subscript[TransMatL, 8])[(Subscript[h, 2])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[e, 3])];B8=(Subscript[TransMatL, 8])[(Subscript[f, 3])];C8=(Subscript[TransMatL, 8])[(Subscript[e, 4])];D8=(Subscript[TransMatL, 8])[(Subscript[f, 4])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[g, 3])];B8=(Subscript[TransMatL, 8])[(Subscript[h, 3])];C8=(Subscript[TransMatL, 8])[(Subscript[g, 4])];D8=(Subscript[TransMatL, 8])[(Subscript[h, 4])];D16=BuildUp[8,A8,B8,C8,D8];BB=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatL, 8])[(Subscript[a, 5])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 5])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 6])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 6])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[c, 5])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 5])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 6])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 6])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[a, 7])];B8=(Subscript[TransMatL, 8])[(Subscript[b, 7])];C8=(Subscript[TransMatL, 8])[(Subscript[a, 8])];D8=(Subscript[TransMatL, 8])[(Subscript[b, 8])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[c, 7])];B8=(Subscript[TransMatL, 8])[(Subscript[d, 7])];C8=(Subscript[TransMatL, 8])[(Subscript[c, 8])];D8=(Subscript[TransMatL, 8])[(Subscript[d, 8])];D16=BuildUp[8,A8,B8,C8,D8];CC=BuildUp[16,A16,B16,C16,D16];

A8=(Subscript[TransMatL, 8])[(Subscript[e, 5])];B8=(Subscript[TransMatL, 8])[(Subscript[f, 5])];C8=(Subscript[TransMatL, 8])[(Subscript[e, 6])];D8=(Subscript[TransMatL, 8])[(Subscript[f, 6])];A16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[g, 5])];B8=(Subscript[TransMatL, 8])[(Subscript[h, 5])];C8=(Subscript[TransMatL, 8])[(Subscript[g, 6])];D8=(Subscript[TransMatL, 8])[(Subscript[h, 6])];B16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[e, 7])];B8=(Subscript[TransMatL, 8])[(Subscript[f, 7])];C8=(Subscript[TransMatL, 8])[(Subscript[e, 8])];D8=(Subscript[TransMatL, 8])[(Subscript[f, 8])];C16=BuildUp[8,A8,B8,C8,D8];A8=(Subscript[TransMatL, 8])[(Subscript[g, 7])];B8=(Subscript[TransMatL, 8])[(Subscript[h, 7])];C8=(Subscript[TransMatL, 8])[(Subscript[g, 8])];D8=(Subscript[TransMatL, 8])[(Subscript[h, 8])];D16=BuildUp[8,A8,B8,C8,D8];DD=BuildUp[16,A16,B16,C16,D16];

st=AbsoluteTime[];(Subscript[PEP, 4])[AA,BB,CC,DD] ;AbsoluteTime[]-st  (*  This takes less than half a dozen seconds.  *)

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[32]&&AA.X+BB.Z==ZeroMat[32]&&CC.W+DD.Y==ZeroMat[32]&&CC.X+DD.Z==IdentityMatrix[32]

{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4])}

Floor[N[W.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+X.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]]]

{-1,-9,4,8,-8,0,-8,6,     -9,2,-9,-3,0,-2,-4,-5,     -3,-7,1,-8,6,-9,5,2,     -8,-3,2,-3,-3,1,2,-7}

Y.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+Z.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]==Join[(Subscript[x, 5]),(Subscript[x, 6]),(Subscript[x, 7]),(Subscript[x, 8])]

(*  With there being 64 coefficients, there exist 2 possibilities for each coefficient. *)        (* Let's see that's about how many distinct systems that have these coefficients? ...  *)
(2^64)

18446744073709551616         
(*  "Oh my," said Dorothy, "Oh my, there are quotro-lions and octo-tigers and hexa-bears."  *)

(*  We've done 3 of them 18-kajillion. Mercifully, I think that's enough. Dorothy agrees.  *)

(*  ------------------------------------------  fin: large octonion systems  -----------------------------------  *)

(*  To illustrate the solution we randomly select 16 sedenion constants with +/- digit components .  *)

SeedRandom[9797];(Subscript[a, 1])=TFR[16];(Subscript[a, 2])=TFR[16];TFR[16];(Subscript[b, 1])=TFR[16];(Subscript[b, 2])=TFR[16];(Subscript[c, 1])=TFR[16];(Subscript[c, 2])=TFR[16];(Subscript[d, 1])=TFR[16];(Subscript[d, 2])=TFR[16];(Subscript[a, 3])=TFR[16];(Subscript[a, 4])=TFR[16];(Subscript[b, 3])=TFR[16];(Subscript[b, 4])=TFR[16];(Subscript[c, 3])=TFR[16];(Subscript[c, 4])=TFR[16];(Subscript[d, 3])=TFR[16];(Subscript[d, 4])=TFR[16];
(*  User must "ENTER" here. You may (of course) change the seed number above to generate new systems.  *)

(*   As is often done in textbook examples, randomly select the answers so that the solution will avoid 
those ugly fractions.  *)

(*  "Answers" for (Subscript[x, 1]), ..., (Subscript[x, 4])  *)(*  User must "ENTER" below.  *)

(Subscript[x, 1])=TFR[16] ;(Subscript[x, 2])=TFR[16] ;(Subscript[x, 3])=TFR[16] ;(Subscript[x, 4])=TFR[16];

(*  Compute the rhs constants for multiplication from the left.  *)

(Subscript[k, 1])=SedDot[(Subscript[a, 1]),(Subscript[x, 1])]+SedDot[(Subscript[b, 1]),(Subscript[x, 2])]+SedDot[(Subscript[c, 1]),(Subscript[x, 3])]+SedDot[(Subscript[d, 1]),(Subscript[x, 4])]           (*  User must "ENTER" here.  *)

(Subscript[k, 2])=SedDot[(Subscript[a, 2]),(Subscript[x, 1])]+SedDot[(Subscript[b, 2]),(Subscript[x, 2])] +SedDot[(Subscript[c, 2]),(Subscript[x, 3])]+SedDot[(Subscript[d, 2]),(Subscript[x, 4])]           (*  User must "ENTER" here.  *)

(Subscript[k, 3])=SedDot[(Subscript[a, 3]),(Subscript[x, 1])]+SedDot[(Subscript[b, 3]),(Subscript[x, 2])]+SedDot[(Subscript[c, 3]),(Subscript[x, 3])]+SedDot[(Subscript[d, 3]),(Subscript[x, 4])]            (*  User must "ENTER" here.  *)

(Subscript[k, 4])=SedDot[(Subscript[a, 4]),(Subscript[x, 1])]+SedDot[(Subscript[b, 4]),(Subscript[x, 2])]+SedDot[(Subscript[c, 4]),(Subscript[x, 3])]+SedDot[(Subscript[d, 4]),(Subscript[x, 4])]            (*  User must "ENTER" here.  *)

(*     We solve the sedenion system:
	 (Subscript[a, 1])(Subscript[x, 1]) + (Subscript[b, 1])(Subscript[x, 2]) + (Subscript[c, 1])(Subscript[x, 3]) + (Subscript[d, 1])(Subscript[x, 4]) = (Subscript[k, 1]),
		(Subscript[a, 2])(Subscript[x, 1]) + (Subscript[b, 2])(Subscript[x, 2]) + (Subscript[c, 2])(Subscript[x, 3]) + (Subscript[d, 2])(Subscript[x, 4]) = (Subscript[k, 2]),
		(Subscript[a, 3])(Subscript[x, 1]) + (Subscript[b, 3])(Subscript[x, 2]) + (Subscript[c, 3])(Subscript[x, 3]) + (Subscript[d, 3])(Subscript[x, 4]) = (Subscript[k, 1]),
		(Subscript[a, 4])(Subscript[x, 1]) + (Subscript[b, 4])(Subscript[x, 2]) + (Subscript[c, 4])(Subscript[x, 3]) + (Subscript[d, 4])(Subscript[x, 4]) = (Subscript[k, 1]),		
  *)

A16=(Subscript[TransMatL, 16])[(Subscript[a, 1])];B16=(Subscript[TransMatL, 16])[(Subscript[b, 1])];C16=(Subscript[TransMatL, 16])[(Subscript[a, 2])];D16=(Subscript[TransMatL, 16])[(Subscript[b, 2])];AA=BuildUp[16,A16,B16,C16,D16];

A16=(Subscript[TransMatL, 16])[(Subscript[c, 1])];B16=(Subscript[TransMatL, 16])[(Subscript[d, 1])];C16=(Subscript[TransMatL, 16])[(Subscript[c, 2])];D16=(Subscript[TransMatL, 16])[(Subscript[d, 2])];BB=BuildUp[16,A16,B16,C16,D16];

A16=(Subscript[TransMatL, 16])[(Subscript[a, 3])];B16=(Subscript[TransMatL, 16])[(Subscript[b, 3])];C16=(Subscript[TransMatL, 16])[(Subscript[a, 4])];D16=(Subscript[TransMatL, 16])[(Subscript[b, 4])];CC=BuildUp[16,A16,B16,C16,D16];

A16=(Subscript[TransMatL, 16])[(Subscript[c, 3])];B16=(Subscript[TransMatL, 16])[(Subscript[d, 3])];C16=(Subscript[TransMatL, 16])[(Subscript[c, 4])];D16=(Subscript[TransMatL, 16])[(Subscript[d, 4])];DD=BuildUp[16,A16,B16,C16,D16];

(Subscript[PEP, 2])[AA,BB,CC,DD]

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[32]&&AA.X+BB.Z==ZeroMat[32]&&CC.W+DD.Y==ZeroMat[32]&&CC.X+DD.Z==IdentityMatrix[32]

{(Subscript[x, 1]),(Subscript[x, 2])}

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]

{0,1,-9,-9,3,-2,6,0,0,5,-8,-6,8,7,1,5,       9,-8,8,-8,7,8,2,-7,-3,-2,2,-7,3,9,5,3}

Y.Join[(Subscript[k, 1]),(Subscript[k, 2])]+Z.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 3]),(Subscript[x, 4])]

(*  With all sedenion coefficients from the right, ... we illustrate solving the octonion system:

	 (Subscript[x, 1])(Subscript[a, 1]) + (Subscript[x, 2])(Subscript[b, 1]) + (Subscript[x, 3])(Subscript[c, 1]) + (Subscript[x, 4])(Subscript[d, 1]) = (Subscript[k, 1]),
		(Subscript[x, 1])(Subscript[a, 2]) + (Subscript[x, 2])(Subscript[b, 2]) + (Subscript[x, 3])(Subscript[c, 2]) + (Subscript[x, 4])(Subscript[d, 2]) = (Subscript[k, 2]),
		(Subscript[x, 1])(Subscript[a, 3]) + (Subscript[x, 2])(Subscript[b, 3]) + (Subscript[x, 3])(Subscript[c, 3]) + (Subscript[x, 4])(Subscript[d, 3]) = (Subscript[k, 3]),
		(Subscript[x, 1])(Subscript[a, 4]) + (Subscript[x, 2])(Subscript[b, 4]) + (Subscript[x, 3])(Subscript[c, 4]) + (Subscript[x, 4])(Subscript[d, 4]) = (Subscript[k, 4]),		
  *)

(Subscript[k, 1])=SedDot[(Subscript[x, 1]),(Subscript[a, 1])]+SedDot[(Subscript[x, 2]),(Subscript[b, 1])]+SedDot[(Subscript[x, 3]),(Subscript[c, 1])]+SedDot[(Subscript[x, 4]),(Subscript[d, 1])]           (*  User must "ENTER" here.  *)

(Subscript[k, 2])=SedDot[(Subscript[x, 1]),(Subscript[a, 2])]+SedDot[(Subscript[x, 2]),(Subscript[b, 2])] +SedDot[(Subscript[x, 3]),(Subscript[c, 2])]+SedDot[(Subscript[x, 4]),(Subscript[d, 2])]           (*  User must "ENTER" here.  *)

(Subscript[k, 3])=SedDot[(Subscript[x, 1]),(Subscript[a, 3])]+SedDot[(Subscript[x, 2]),(Subscript[b, 3])]+SedDot[(Subscript[x, 3]),(Subscript[c, 3])]+SedDot[(Subscript[x, 4]),(Subscript[d, 3])]            (*  User must "ENTER" here.  *)

(Subscript[k, 4])=SedDot[(Subscript[x, 1]),(Subscript[a, 4])]+SedDot[(Subscript[x, 2]),(Subscript[b, 4])]+SedDot[(Subscript[x, 3]),(Subscript[c, 4])]+SedDot[(Subscript[x, 4]),(Subscript[d, 4])]            (*  User must "ENTER" here.  *)

A16=(Subscript[TransMatR, 16])[(Subscript[a, 1])];B16=(Subscript[TransMatR, 16])[(Subscript[b, 1])];C16=(Subscript[TransMatR, 16])[(Subscript[a, 2])];D16=(Subscript[TransMatR, 16])[(Subscript[b, 2])];AA=BuildUp[16,A16,B16,C16,D16];

A16=(Subscript[TransMatR, 16])[(Subscript[c, 1])];B16=(Subscript[TransMatR, 16])[(Subscript[d, 1])];C16=(Subscript[TransMatR, 16])[(Subscript[c, 2])];D16=(Subscript[TransMatR, 16])[(Subscript[d, 2])];BB=BuildUp[16,A16,B16,C16,D16];

A16=(Subscript[TransMatR, 16])[(Subscript[a, 3])];B16=(Subscript[TransMatR, 16])[(Subscript[b, 3])];C16=(Subscript[TransMatR, 16])[(Subscript[a, 4])];D16=(Subscript[TransMatR, 16])[(Subscript[b, 4])];
CC=BuildUp[16,A16,B16,C16,D16];

A16=(Subscript[TransMatR, 16])[(Subscript[c, 3])];B16=(Subscript[TransMatR, 16])[(Subscript[d, 3])];C16=(Subscript[TransMatR, 16])[(Subscript[c, 4])];D16=(Subscript[TransMatR, 16])[(Subscript[d, 4])];DD=BuildUp[16,A16,B16,C16,D16];

(Subscript[PEP, 4])[AA,BB,CC,DD]

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[32]&&AA.X+BB.Z==ZeroMat[32]&&CC.W+DD.Y==ZeroMat[32]&&CC.X+DD.Z==IdentityMatrix[32]

{(Subscript[x, 1]),(Subscript[x, 2])}

W.Join[(Subscript[k, 1]),(Subscript[k, 2])]+X.Join[(Subscript[k, 3]),(Subscript[k, 4])]

{0,1,-9,-9,3,-2,6,0,0,5,-8,-6,8,7,1,5,        9,-8,8,-8,7,8,2,-7,-3,-2,2,-7,3,9,5,3}

Y.Join[(Subscript[k, 1]),(Subscript[k, 2])]+Z.Join[(Subscript[k, 3]),(Subscript[k, 4])]==Join[(Subscript[x, 3]),(Subscript[x, 4])]

(*  We select 64 constants as coefficients of an 8x8 system of sedenions.  *)

SeedRandom[9797];h=16;(Subscript[a, 1])=TFR[h];(Subscript[a, 2])=TFR[h];(Subscript[b, 1])=TFR[h];(Subscript[b, 2])=TFR[h];(Subscript[c, 1])=TFR[h];(Subscript[c, 2])=TFR[h];(Subscript[d, 1])=TFR[h];(Subscript[d, 2])=TFR[h];(Subscript[a, 3])=TFR[h];(Subscript[a, 4])=TFR[h];(Subscript[b, 3])=TFR[h];(Subscript[b, 4])=TFR[h];(Subscript[c, 3])=TFR[h];(Subscript[c, 4])=TFR[h];(Subscript[d, 3])=TFR[h];(Subscript[d, 4])=TFR[h];(Subscript[e, 1])=TFR[h];(Subscript[f, 1])=TFR[h];(Subscript[g, 1])=TFR[h];(Subscript[h, 1])=TFR[h];(Subscript[e, 2])=TFR[h];(Subscript[f, 2])=TFR[h];(Subscript[g, 2])=TFR[h];(Subscript[h, 2])=TFR[h];(Subscript[e, 3])=TFR[h];(Subscript[f, 3])=TFR[h];(Subscript[g, 3])=TFR[h];(Subscript[h, 3])=TFR[h];(Subscript[e, 4])=TFR[h];(Subscript[f, 4])=TFR[h];(Subscript[g, 4])=TFR[h];(Subscript[h, 4])=TFR[h];(Subscript[e, 5])=TFR[h];(Subscript[f, 5])=TFR[h];(Subscript[g, 5])=TFR[h];(Subscript[h, 5])=TFR[h];(Subscript[e, 6])=TFR[h];(Subscript[f, 6])=TFR[h];(Subscript[g, 6])=TFR[h];(Subscript[h, 6])=TFR[h];(Subscript[e, 7])=TFR[h];(Subscript[f, 7])=TFR[h];(Subscript[g, 7])=TFR[h];(Subscript[h, 7])=TFR[h];(Subscript[e, 8])=TFR[h];(Subscript[f, 8])=TFR[h];(Subscript[g, 8])=TFR[h];(Subscript[h, 8])=TFR[h];(Subscript[a, 5])=TFR[h];(Subscript[b, 5])=TFR[h];(Subscript[c, 5])=TFR[h];(Subscript[d, 5])=TFR[h];(Subscript[a, 6])=TFR[h];(Subscript[b, 6])=TFR[h];(Subscript[c, 6])=TFR[h];(Subscript[d, 6])=TFR[h];(Subscript[a, 7])=TFR[h];(Subscript[b, 7])=TFR[h];(Subscript[c, 7])=TFR[h];(Subscript[d, 7])=TFR[h];(Subscript[a, 8])=TFR[h];(Subscript[b, 8])=TFR[h];(Subscript[c, 8])=TFR[h];(Subscript[d, 8])=TFR[h];
(*  User must "ENTER" here. You may (of course) change the seed number above to generate new systems.  *)

(Subscript[x, 8])=TFR[h];(Subscript[x, 7])=TFR[h];(Subscript[x, 6])=TFR[h];(Subscript[x, 5])=TFR[h];(Subscript[x, 4])=TFR[h];(Subscript[x, 3])=TFR[h];(Subscript[x, 2])=TFR[h];(Subscript[x, 1])=TFR[h];

(*  We have the system with constant multiplying from the right side.  *)

(*  User must "ENTER" here.  *) (Subscript[k, 1])=SedDot[(Subscript[x, 1]),(Subscript[a, 1])]+SedDot[(Subscript[x, 2]),(Subscript[b, 1])]+SedDot[(Subscript[x, 3]),(Subscript[c, 1])]+SedDot[(Subscript[x, 4]),(Subscript[d, 1])]+SedDot[(Subscript[x, 5]),(Subscript[e, 1])]+SedDot[(Subscript[x, 6]),(Subscript[f, 1])]+SedDot[(Subscript[x, 7]),(Subscript[g, 1])]+SedDot[(Subscript[x, 8]),(Subscript[h, 1])];     
(Subscript[k, 2])=SedDot[(Subscript[x, 1]),(Subscript[a, 2])]+SedDot[(Subscript[x, 2]),(Subscript[b, 2])] +SedDot[(Subscript[x, 3]),(Subscript[c, 2])]+SedDot[(Subscript[x, 4]),(Subscript[d, 2])]+SedDot[(Subscript[x, 5]),(Subscript[e, 2])]+SedDot[(Subscript[x, 6]),(Subscript[f, 2])]+SedDot[(Subscript[x, 7]),(Subscript[g, 2])]+SedDot[(Subscript[x, 8]),(Subscript[h, 2])];
(Subscript[k, 3])=SedDot[(Subscript[x, 1]),(Subscript[a, 3])]+SedDot[(Subscript[x, 2]),(Subscript[b, 3])]+SedDot[(Subscript[x, 3]),(Subscript[c, 3])]+SedDot[(Subscript[x, 4]),(Subscript[d, 3])]+SedDot[(Subscript[x, 5]),(Subscript[e, 3])]+SedDot[(Subscript[x, 6]),(Subscript[f, 3])]+SedDot[(Subscript[x, 7]),(Subscript[g, 3])]+SedDot[(Subscript[x, 8]),(Subscript[h, 3])];
(Subscript[k, 4])=SedDot[(Subscript[x, 1]),(Subscript[a, 4])]+SedDot[(Subscript[x, 2]),(Subscript[b, 4])]+SedDot[(Subscript[x, 3]),(Subscript[c, 4])]+SedDot[(Subscript[x, 4]),(Subscript[d, 4])]+SedDot[(Subscript[x, 5]),(Subscript[e, 4])]+SedDot[(Subscript[x, 6]),(Subscript[f, 4])]+SedDot[(Subscript[x, 7]),(Subscript[g, 4])]+SedDot[(Subscript[x, 8]),(Subscript[h, 4])];
(Subscript[k, 5])=SedDot[(Subscript[x, 1]),(Subscript[a, 5])]+SedDot[(Subscript[x, 2]),(Subscript[b, 5])]+SedDot[(Subscript[x, 3]),(Subscript[c, 5])]+SedDot[(Subscript[x, 4]),(Subscript[d, 5])]+SedDot[(Subscript[x, 5]),(Subscript[e, 5])]+SedDot[(Subscript[x, 6]),(Subscript[f, 5])]+SedDot[(Subscript[x, 7]),(Subscript[g, 5])]+SedDot[(Subscript[x, 8]),(Subscript[h, 5])];
(Subscript[k, 6])=SedDot[(Subscript[x, 1]),(Subscript[a, 6])]+SedDot[(Subscript[x, 2]),(Subscript[b, 6])] +SedDot[(Subscript[x, 3]),(Subscript[c, 6])]+SedDot[(Subscript[x, 4]),(Subscript[d, 6])]+SedDot[(Subscript[x, 5]),(Subscript[e, 6])]+SedDot[(Subscript[x, 6]),(Subscript[f, 6])] +SedDot[(Subscript[x, 7]),(Subscript[g, 6])]+SedDot[(Subscript[x, 8]),(Subscript[h, 6])];
(Subscript[k, 7])=SedDot[(Subscript[x, 1]),(Subscript[a, 7])]+SedDot[(Subscript[x, 2]),(Subscript[b, 7])]+SedDot[(Subscript[x, 3]),(Subscript[c, 7])]+SedDot[(Subscript[x, 4]),(Subscript[d, 7])]+SedDot[(Subscript[x, 5]),(Subscript[e, 7])]+SedDot[(Subscript[x, 6]),(Subscript[f, 7])]+SedDot[(Subscript[x, 7]),(Subscript[g, 7])]+SedDot[(Subscript[x, 8]),(Subscript[h, 7])];
(Subscript[k, 8])=SedDot[(Subscript[x, 1]),(Subscript[a, 8])]+SedDot[(Subscript[x, 2]),(Subscript[b, 8])]+SedDot[(Subscript[x, 3]),(Subscript[c, 8])]+SedDot[(Subscript[x, 4]),(Subscript[d, 8])]+SedDot[(Subscript[x, 5]),(Subscript[e, 8])]+SedDot[(Subscript[x, 6]),(Subscript[f, 8])]+SedDot[(Subscript[x, 7]),(Subscript[g, 8])]+SedDot[(Subscript[x, 8]),(Subscript[h, 8])];

{(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4]),(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])}

Ah=(Subscript[TransMatR, 16])[(Subscript[a, 1])];Bh=(Subscript[TransMatR, 16])[(Subscript[b, 1])];Ch=(Subscript[TransMatR, 16])[(Subscript[a, 2])];Dh=(Subscript[TransMatR, 16])[(Subscript[b, 2])];A32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[c, 1])];Bh=(Subscript[TransMatR, 16])[(Subscript[d, 1])];Ch=(Subscript[TransMatR, 16])[(Subscript[c, 2])];Dh=(Subscript[TransMatR, 16])[(Subscript[d, 2])];B32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[a, 3])];Bh=(Subscript[TransMatR, 16])[(Subscript[b, 3])];Ch=(Subscript[TransMatR, 16])[(Subscript[a, 4])];Dh=(Subscript[TransMatR, 16])[(Subscript[b, 4])];C32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[c, 3])];Bh=(Subscript[TransMatR, 16])[(Subscript[d, 3])];Ch=(Subscript[TransMatR, 16])[(Subscript[c, 4])];Dh=(Subscript[TransMatR, 16])[(Subscript[d, 4])];D32=BuildUp[h,Ah,Bh,Ch,Dh];AA=BuildUp[32,A32,B32,C32,D32];

Ah=(Subscript[TransMatR, 16])[(Subscript[e, 1])];Bh=(Subscript[TransMatR, 16])[(Subscript[f, 1])];Ch=(Subscript[TransMatR, 16])[(Subscript[e, 2])];Dh=(Subscript[TransMatR, 16])[(Subscript[f, 2])];A32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[g, 1])];Bh=(Subscript[TransMatR, 16])[(Subscript[h, 1])];Ch=(Subscript[TransMatR, 16])[(Subscript[g, 2])];Dh=(Subscript[TransMatR, 16])[(Subscript[h, 2])];B32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[e, 3])];Bh=(Subscript[TransMatR, 16])[(Subscript[f, 3])];Ch=(Subscript[TransMatR, 16])[(Subscript[e, 4])];Dh=(Subscript[TransMatR, 16])[(Subscript[f, 4])];C32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[g, 3])];Bh=(Subscript[TransMatR, 16])[(Subscript[h, 3])];Ch=(Subscript[TransMatR, 16])[(Subscript[g, 4])];Dh=(Subscript[TransMatR, 16])[(Subscript[h, 4])];D32=BuildUp[h,Ah,Bh,Ch,Dh];BB=BuildUp[32,A32,B32,C32,D32];

Ah=(Subscript[TransMatR, 16])[(Subscript[a, 5])];Bh=(Subscript[TransMatR, 16])[(Subscript[b, 5])];Ch=(Subscript[TransMatR, 16])[(Subscript[a, 6])];Dh=(Subscript[TransMatR, 16])[(Subscript[b, 6])];A32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[c, 5])];Bh=(Subscript[TransMatR, 16])[(Subscript[d, 5])];Ch=(Subscript[TransMatR, 16])[(Subscript[c, 6])];Dh=(Subscript[TransMatR, 16])[(Subscript[d, 6])];B32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[a, 7])];Bh=(Subscript[TransMatR, 16])[(Subscript[b, 7])];Ch=(Subscript[TransMatR, 16])[(Subscript[a, 8])];Dh=(Subscript[TransMatR, 16])[(Subscript[b, 8])];C32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[c, 7])];Bh=(Subscript[TransMatR, 16])[(Subscript[d, 7])];Ch=(Subscript[TransMatR, 16])[(Subscript[c, 8])];Dh=(Subscript[TransMatR, 16])[(Subscript[d, 8])];D32=BuildUp[h,Ah,Bh,Ch,Dh];CC=BuildUp[32,A32,B32,C32,D32];

Ah=(Subscript[TransMatR, 16])[(Subscript[e, 5])];Bh=(Subscript[TransMatR, 16])[(Subscript[f, 5])];Ch=(Subscript[TransMatR, 16])[(Subscript[e, 6])];Dh=(Subscript[TransMatR, 16])[(Subscript[f, 6])];A32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[g, 5])];Bh=(Subscript[TransMatR, 16])[(Subscript[h, 5])];Ch=(Subscript[TransMatR, 16])[(Subscript[g, 6])];Dh=(Subscript[TransMatR, 16])[(Subscript[h, 6])];B32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[e, 7])];Bh=(Subscript[TransMatR, 16])[(Subscript[f, 7])];Ch=(Subscript[TransMatR, 16])[(Subscript[e, 8])];Dh=(Subscript[TransMatR, 16])[(Subscript[f, 8])];C32=BuildUp[h,Ah,Bh,Ch,Dh];Ah=(Subscript[TransMatR, 16])[(Subscript[g, 7])];Bh=(Subscript[TransMatR, 16])[(Subscript[h, 7])];Ch=(Subscript[TransMatR, 16])[(Subscript[g, 8])];Dh=(Subscript[TransMatR, 16])[(Subscript[h, 8])];D32=BuildUp[h,Ah,Bh,Ch,Dh];DD=BuildUp[32,A32,B32,C32,D32];

st=AbsoluteTime[];(Subscript[PEP, 1])[AA,BB,CC,DD] ;AbsoluteTime[]-st  (*  This takes less than a minute and a half.  *)

(*  Showing that:        ({{A, B}, {C, D}}) * ({{W, X}, {Y, Z}})   =   ({{I, O}, {O, I}})              *)

AA.W+BB.Y==IdentityMatrix[64]&&AA.X+BB.Z==ZeroMat[64]&&CC.W+DD.Y==ZeroMat[64]&&CC.X+DD.Z==IdentityMatrix[64]

Join[(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4])]

W.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+X.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]

Y.Join[(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3]),(Subscript[k, 4])]+Z.Join[(Subscript[k, 5]),(Subscript[k, 6]),(Subscript[k, 7]),(Subscript[k, 8])]==Join[(Subscript[x, 5]),(Subscript[x, 6]),(Subscript[x, 7]),(Subscript[x, 8])]

(*  --------------------------------------------------  FIN  -----------------------------------------------------------  *)